I. Introduction

Since the pioneering work of Bernoulli [1], mathematical models have proved to be a fundamental tool for predicting the dynamics of infectious diseases. In particular, the recent changes in the UK’s strategies regarding the COVID-19 pandemic [2] strongly collaborate to reaffirm such a claim. It turns out that an enormous number of mathematical models and simulating scenarios were recently proposed and analyzed for the COVID-19 outbreak, making it hard to attempt a complete literature overview. Nevertheless, at the heart of many of these proposed models are SIR-like compartmental models, originally proposed by Kermack & Mckendrick [3]. In other words, unless some compartment clustering (summarizing individuals in distinct compartments) is used, such proposed models assume that at any time $t\ge 0$ of the underlying diseases, the total population $\textcolor{red}{}N\left(t\right)=S_{\mathrm{(t)}}+I_{\mathrm{(t)}}+R_{(t)}$is divided into compartments of susceptible $S_{\mathrm{(t)}}$,infected or infectious $I_{\mathrm{(t)}}$and recovered or removed $R_{\mathrm{(t)}}$.

Using the law of mass action, Kermack & Mckendrick [3] concludes that the diseases must follow the SIR nonlinear dynamics:

(1)

where $\text{'}:=\frac{d}{dt},\beta$is a disease transmission rate, $\mu$ is a birth and mortality rates (that for simplicity we assume to be the same), $\gamma$ is the inverse of the mean infectious period. We assume that all parameters are time-independent.

In the dynamics (1), the following initial conditions are considered (IC)

$S\mathrm{(0)}=1-I(0),I(0)\ge 0,R(0)=\; 0,(2)$

characterizing the model as initial value problems (IVP) for nonlinear systems of ordinary differential equations (ODE’s).

Given the non-linearity, SIR model (1) - (2) does not have an explicit solution. As a result, any attempt to predict disease dynamics must be obtained from some approximation of the solution provided by numerical methods [4]. This necessity of using numerical methods drives the main question of this contribution: Are there some parameter choice scenarios that characterizes the SIR model to be stiff?

The phenomenon of stiffness in a (IVP) means that the numerical solutions are unstable and/or its accuracy is largely affected using any explicit numerical methods. It imposes the need for a very small time step to approximate the solution and, consequently, becomes computationally expensive.

In Section II, we revisit some well-known properties of the solution of the (IVP) of the SIR model. Such properties will help to justify the numerical analysis that follows. In Section III, we present a definition of stiffness. Moreover, we present and compare the obtained results for distinct scenarios of parameter choices for the model. In Section IV, we show the consequences of the stiffness phenomena in the numerical solution of the SIR model using explicit and implicit numerical approximations. The main contribution of this manuscript is the conclusion that the SIR model is indeed stiff, for large transmission rates $\beta$. The analysis of the error in the numericla solutions allows us to recommend the use of appropriate numerical methods as well as the time-step size. In Section V, we formulate the final conclusions.

II. Some sir model properties

In this section, we revisit some of the well-known SIR model properties reformulated from [5], that will be used to justify the numerical analysis that follows.

Theorem 1. Let the assumptions for the SIR model (1) with initial conditions (2) hold true. Then, the (IVP) (1) - (2) corresponding to the SIR model have a unique differentiable solution$Y(t)={\mathrm{(S(t),I(t),R(t))}}^t\textcolor{red}{}$

Since the initial conditions (2) for are non-negative, each of the coordinates of the solution $\mathrm{Y(t)}\textcolor{red}{}$ remains non-negative for all $t\; \ge \; 0$

Note that, summing all the equations in models (1), we conclude that the total population is constant. To keep things simple, we will assume that Importar imagen is normalized so that Importar imagen , for all

It follows from the first equation in that is decreasing. Furthermore, we get from Theorem 1 that

(3)

As we will see in Section IV, the stiffness phenomena present in the SIR model imply that numerical solutions obtained from an explicit numerical scheme might not satisfy Theorem 1 and equation (3) in particular.

Note that, summing all the equations in models (1), we conclude that the total population is constant. To keep things simple, we will assume that $N$ is normalized so that $\mathrm{N(t)}=\; 1$ , for all $t\ge 0$.

It follows from the first equation in that Importar imagen is decreasing. Furthermore, we get from Theorem 1 that

(3)

As we will see in Section IV, the stiffness phenomena present in the SIR model imply that numerical solutions obtained from an explicit numerical scheme might not satisfy Theorem 1 and equation (3) in particular.

III. Numerical Methods And Stiffness Phenomena For Sir Model

Let $\mathrm{Y(t)}$ be the unique solution of the SIR model (1) - (2) as shown in Theorem 1, whose coordinates will be denoted by $Y^i\mathrm{(t)}$ for $i=1,2,3$respectively. According to [4], the step of a numerical method to obtain the approximation $Y_n^i$ for the exact solution $Y^i\left(t_n\right)$ of the (IVP), for $t_n=\in [0,T]$, can be written, in general, as follows

where $k\; \ge \; 1$ is the number of preview steps used to compute $Y_n^i$. The values ${\alpha }_j$ , for $\textcolor{red}{}j=0,..,k$ are numerical method-specific constants, Importar imagen is the method step-size. The function Importar imagen is defined in terms of evaluating the right-hand side Importar imagen of the system of differential equations (1).

In general terms, we can classify the numerical methods as: Single-step if $\mathrm{k=1}$, or multi-step if $k>1$ . Linear if ${\psi }_F$ is linear in the evaluation of $\mathrm{F(t,Y(t))}$ and non-linear otherwise. Implicit if ${\psi }_f$ depends on $Y_n$ and explicit otherwise.

In the numerical simulations for the SIR model (1) - (2) that follows, we will adopt the following methods: A) the explicit, single-step Euler method; B) the explicit, single-step Runge-Kutta method; C) the semi-implicit, fourth order, four-stage Rosenbrock method, e.g, [4].

A. The stiffness phenomena

The phenomena of stiffness in solving (IVP) appeared at the beginning of the 20th century. It was [6] one of the first to notice some implicit numerical methods perform better than explicit ones for stiff (IVP). This is because the implicit methods have an unlimited stability region that covers the entire half-plane complex with a negative real part or at least an unlimited part of it, e.g. [4].

A well-accepted definition of the phenomena of stiffness follows:

Definition 1. Let ${\lambda }_j\left(t\right)\in ℂ$, for $j=1,...,n$ the eigenvalues of the Jacobian matrix $\mathrm{J(F(t,Y(t)))}$ , where $\textcolor{red}{}F(t,Y(t\left)\right)$ is the right-hand side of the (IVP).

The (IVP) is called stiff in some interval I of the real line if $\textcolor{red}{}\mathrm{Re(}{\lambda }_j)<0$ for all $\mathrm{j=1,...,n}$ and stiffness ratio

If the (IVP) is stiff, then certain numerical methods for solving the equation become numerically unstable unless the step size is taken to be extremely small. This is the case for most explicit numerical schemes. See [4] for some examples.

B. The stiffness phenomena for SIR model

The Jacobian matrix for model (1) for an arbitrary point is given by

Then, the eigenvalues ${\lambda }_j\mathrm{,\; j=1,2,3}$ of the Jacobian matrix are given by

Below, the stiffness ratio for the (IVP) SIR model (1) - (2) is numerically investigated. Table I displays various values of the stiffness ratio $\xi$ (see Definition 1) as a function of $\beta$ and$\mu$. The numerical results presented in Table I suggest that the SIR model’s rigidity becomes increasingly greater as we increase the $\beta$ and/or decrease the $\mu$.

TABLE I.

Stiffness ratio ξ for SIR model, as a function of $\beta$ and $\mu$

β	μ=0.01	μ=0.001	μ=0.0001
0.5	ξ=12	ξ=131	ξ=1354
1.5	ξ=39	ξ=392	ξ=3923
10	ξ=808	ξ=8090	ξ=80908
12	ξ=986	ξ=9870	ξ=98708
18	ξ=1442	ξ=14430	ξ=144308
25	ξ=1766	ξ=17666	ξ=176668

Before delving into the effect of the stiffness ratio on numerical simulations of the SIR model (1), we should consider the epidemiological implications of the parameters $\mu$ and $\beta$ as shown in Table I.

Remark 1. Small values of $\mu$ mean that the birth and mortality rates are very low, in proportion to the total population, during the disease duration. It is a typical effect for a short time simulation. Large values of the transmission rate $\beta$ have a more debatable epidemiological meaning. However, there are highly contagious diseases such as measles, pertussis, and tuberculosis that are already documented in the literature [7, 8]. Furthermore, many diseases transmission depends on population density (high for high-density populations) and environment (for example, for respiratory transmitted diseases, a crowded and closed environment implies higher transmission rates) [9]. To put it another way, a disease that has not yet been reported may have a higher transmission rate than those already reported.

IV. Numerical Results

In this section, we present numerical simulated results and show the implications of the stiffness ratio (see the stiff ration in Table I on the numerical approximated solutions for the SIR model. We present the simulations using e explicit Euler method, the two-stage and second-order Runge-Kutta method (RK22), the four-stage and fourth-order Runge-Kutta method (RK44) and Rosenbrock method [4]. In the simulation, the step-size $\mathrm{h=0.1}$ is used for all the methods except for the Rosenbrock method, for which the adaptive step-size is used [4]. The inicial conditions $S_o\mathrm{=\; 0.9}$, ${\mathrm{I}}_0\mathrm{=\; 0.1}$ and $\textcolor{red}{}{R}_o=0$ are used and the parameter $\gamma \mathrm{=0.35}$ is keeped fixed. The above mentioned methods are implemented in Fortran 90 with a double precision, GNU Fortran Compiler (Ubuntu 18.04.5 LTS operational systema 64 bits) and running in a processor Intel^R Core (TM)i5-5200U CPU@2.20GHzx4

Fig. 1. Simulated scenarios for the infected populations I(t)
(a) Simulated scenarios using the Euler method.

(b)
Simulated scenarios using the RK22 method.

(c)
Simulated scenarios using the RK44 method.

(d)
Simulated scenarios using the Rosenbrock method

In Fig. 1, we show the simulated scenarios of the infected population, for choices of $\beta$ and values of $\beta$ = 0.5, 1.5, 5.0, 10, 12, 18, 25, respectively. In the simulation with the Rosenbrock method, the adaptive step-size, typically of the method, is used. The simulations for the remaining methods use the step-size h=0.1.

For the simulated scenario with $\mathrm{\beta \; =\; 25}$, the first and most obvious consequences of the stiffness phenomenon can be seen. In this particular scenario, the approximate solution using the Euler method oscillates as shown in Fig. 1(a). The oscillations disappear in the numerical solutions using the RK22, RK44 and Rosenbrock methods, as can be observed in Fig. 1(b), Fig. 1(c) and Fig. 1(d).

Fig. 2.
Simulated scenarios using the numerical methods for B=25.
(a) Susceptible population S(t).

(b) Infected population I(t)

A zoom on the numerically simulated scenario for $\mathrm{\beta \; =\; 25}$of the susceptible and infected is depicted in Fig. 2 (a) and Fig. 2(b), respectively, using all the numerical methods under investigation. It turns out that the oscillations of the numerical solutions obtained by the explicit Euler method produce negative values for the susceptible (see Fig. 2 (a)) and values larger than 1 for the infected population (see Fig. 2 (b)), which contradicts Theorem 1. Moreover, the simulated scenarios using the RK22 and RK44 methods present non-smooth solutions and hence also contradict Theorem 1. It can be observed in the simulated results that only the Rosenbrock method, the method designed for stiff problems [4], presents an approximate solution that satisfies Theorem 1.

The remaining question addressed in this contribution is whether the stiffness phenomenon for the SIR model appears only for very large $\beta$

The simulated solution for the susceptible and infected population, with the choice of $\mathrm{\beta \; =\; 12}$ ,is depicted in Fig.3(a) and Fig. 3(b, using the Euler method with step-size $\mathrm{h\; =\; 0.1}$(see Fig. 3(a)).

Fig. 3. Simulated scenarios using the numerical methods of Euler and Rosenbrock for B=12.
(a) Susceptible population.

(b) Infected population.

In Fig. 4(a), we show the values for the step-sizes used by the Rosenbrock method with adaptive step for the simulated scenarios with$\mathrm{\beta \; =\; 12}$. It is worth mentioning that a small step-size is needed in the regions close to the peak of infection; that coincides with the region of the large discrepancies between the two methods presented in (see Fig. 3). In Fig. 4(b), we depicted the numerical solutions for the susceptible and infected populations using the Euler and Rosenbrock methods with the smallest step-size $h\mathrm{=\; 0.0027}$ adopted by the Rosenbrock method. As a result of using a very small step-size, the Euler method shows a more accurate numerical solution. Indeed, it is a consequence of the stiffness of the SIR model.

Fig. 4. h variable (a); Simulated solutions of S(t) and I(t) using the Euler method with h = 0.0027 and Rosenbrock method (b).
(a) Adaptive step-size h used by the Rosenbrock method.

(b)
(b) Numerical solutions obtained for of S(t) and I(t) by the Euler and Rosenbrock methods using the step-size h = 0.0027 for β = 12.

Table II and Table III show the difference in the maximum and quadratic norm of the simulated infected population using the methods discussed in this contribution for different scenarios of $\beta$ .

The difference between the numerical solutions of the infected population using the Euler and Rosenbrock methods with step-size $h\mathrm{=\; 0.1}$is depicted in Fig. 5 for some scenarios of choice for $\beta$ . The absolute value of the difference is monotonically increases with $\beta$, as shown in the second column of Table II. The difference is approximately 10% n in absolute value for $\mathrm{\beta \; =\; 10}$, 20% in absolute values for $\mathrm{\beta \; =\; 18}$, and 40% in absolute values for $\mathrm{\beta \; =25}$.

The difference between the numerical solutions of the infected population using the RK44 and Rosenbrock methods with step-size $h\mathrm{=\; 0.1}$ is depicted in Fig. 6 for some scenarios of choices of $\beta$. As it can be seen from the results in the fourth column of Table II, the difference is significant only for $\mathrm{\beta \; =\; 25}$(see also Fig. 5(d)).

Based on the results in Fig. 5, Fig. 6 and the one in the Table II and Table III, we can conclude that as the stiffness of the SIR model increases (as $\beta$ increases), the absolute difference between the Euler or RK44 and Rosenbrock methods becomes concentrated near the infection’s peak. The results show that the stiffness of the SIR model suggests using a small numerical step-size or a method designed for solving stiff problems, such as the Rosenbrock method.

TABLE II.
ABSOLUTE DIFFERENCE IN THE MAXIMUM NORM E, OBTAINED BY COMPARE THE NUMERICAL SOLUTION OF THE INFECTED POPULATIONS USING EULER, RK22, RK44 AND ROSENBROCK WITH STEP-SIZE h=0.1 FOR DISTINCT SCENARIOS OF CHOICES OF B

, WHERE $E_{\infty }={\max }_{1\le i\le n_{\text{p}}}\left(\right|e_i\left|\right),e_i$ IS THE DIFFERENCE BETWEEN THE NUMERICAL SOLUTION OBTAINED WITH TWO METHODS AND $n_p$ IS THE NUMBER OF MESH POINTS.

β	EUL-RK44	ROSEN-EUL	ROSEN-RK22	ROSEN-RK44
0.5	6.31629999E-004	6.32139999E-004	3.10000000E-006	8.40000000E-007
1.5	9.31865999E-003	9.37654999E-003	6.65800000E-005	8.13800000E-005
5.0	4.31432800E-002	4.45318200E-002	3.14857000E-003	1.53705999E-003
10	9.75375600E-002	0.1046652200000	1.60849599E-002	7.12765999E-003
12	0.1254940800000	0.1386138800000	2.24853200E-002	1.31198000E-002
18	0.1652883499999	0.1921446000000	6.42736900E-002	3.78952900E-002
25	0.25355632000	0.3798633000000	0.1938273099999	0.16452352999999

Fig. 6.
Difference between the numerical solution of the infected population using RK44 and Rosenbrock methods with step-size h = 0.1. Simulated scenarios with distinct choices of β.

(a) β = 0.5

(b)
(b) β = 1.5

(c)
(c) β = 18

(d) β = 25
(d) β = 25

V. CONCLUSIONS

In this work, we numerically investigate the stiffness phenomena for the epidemiological SIR model. The simulated scenarios presented show that the stiffness ratio of the SIR model becomes more pronounced as the disease contagion rate increases or the rate of birth/mortality approaches zero. However, the numerical findings indicate that the explicit numerical method have more difficulties to obtain a numerical solution satisfying Theorem 1 as $\beta$ increases. In particular, the stiffness phenomena imply that the numerical solution using explicit methods, such as Euler and the RK, does not agree with the theoretical results on Theorem 1, while the semi-implicit Rosenbrock method does agree with the theoretical results. The analysis of the maximum and quadratic difference between the analyzed explicit methods and the Rosenbrock method becomes large near the peak of the infection. As a result, the explicit methods overestimate the peak of infection.

We also demonstrate how to use the Rosenbrock method to determine an optimal step-size in terms of $\beta$ to use in the explicit numerical method, resulting in a more accurate numerical solution for the SIR model. This step-size can be used in explicit methods, e.g., the Euler’s method, which is a good alternative to reduce the computational cost.

The results obtained in this contribution collaborate with the following conjecture regarding the use of explicit numerical methods for the SIR model: For a given $\beta$ choose $h$ such that $\beta h<\mathrm{0.1}$ . It will be investigated in detail in future contributions.

References

[1] D. Bernoulli, “Essai D’une Nouvelle Analyse de la Mortalité Causée par la Petite Vérole, et Désavantages de L’inoculation Pour la Prévenir, Histoire de l’Acad., Roy. Sci., Paris, pp. 1-45, 1760.

[2] N. M. Ferguson, D. Laydon, G. Nedjati-Gilani and et al., “Impact of Non-pharmaceutical Interventions (NPIs) to Reduce COVID-19 Mortality and Healthcare Demand”, Imperial College COVID-19 Response Team, London, pp. 1-20, 2020.

[3] W.O. Kermack, and A. G. Mckendrick, “A Contribution to the Mathematical Theory of Epidemics”, Proceedings of the Royal Society of London, London, vol. 115, pp. 700-721, 1927.

[4] J. D. Lambert, “Numerical Methods for Ordinary Differential Systems: The Initial Value Problem”, 1st ed., New York: Wiley, 1991.

[5] H. W. Hethcote, “The Mathematics of Infectious Diseases”, SIAM review, vol. 42, pp. 599-653, 2000.

[6] C. F. Curtiss, and J. O. Hirschfelder, “Integration of Stiff Equations”, Proceedings of the National Academy of Sciences of the United States of America, vol. 38, pp. 235–243, 1952.

[7] J. Kumate, “Infectious Diseases in the 21st Century”, Archives of Medical Research, vol. 28, pp. 155-161, 1997.

[8] I. F. Souza, “Queda nas Taxas de Cobertura Vacinal e Recrudescimento do Sarampo no Brasil: Revisão de Literatura. 2021”, Anais. Bauru: Faculdade de Odontologia de Bauru, Universidade de São Paulo. Available: https://repositorio.usp.br/directbitstream/41cdd6d5-ca58-48b8-9b35-c5c2000c447b/3083766.pdf [Accessed: Aug. 05, 2022], 2021.

[9] J. C. Marques, A. De Cezaro, and M. J. Lazo, “A SIR Model with Spatially Distributed Multiple Populations Interactions for Disease Dissemination”. Trends in Computational and Applied Mathematics - SBMAC, vol. 23, pp. 143-154, 2022.