Estimation of COVID-19 dynamics "on a back-of-envelope": Does the simplest SIR model provide quantitative parameters and predictions?

Basing on existence of the mathematically sequential reduction of the three-compartmental (Susceptible-Infected-Recovered/Removed) model to the Verhulst (logistic) equation with the parameters determined by the basic characteristic of epidemic process, this model is tested in application to the recent data on COVID-19 outbreak reported by the European Centre for Disease Prevention and Control. It is shown that such a simple model adequately reproduces the epidemic dynamics not only qualitatively but for a number of countries quantitatively with a high degree of correlation that allows to use it for predictive estimations. In addition, some features of SIR model are discussed in the context, how its parameters and conditions reflect measures attempted for the disease growth prevention that is also clearly indicated by deviations from such model solutions.

Introduction

The outbreak of disease named by the World Health Organisation (WHO) as COVID-19 induced by the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is now an officially declared as pandemics, which affected practically all world countries, see [1], [2], [3] for a general review.

This situation induces also a demand to the mathematical epidemiologist community [4] for revealing models of outbreak dynamics, which have not only explanatory but also a predictive potential while the outbreak is in an active phase. They are aimed at the fast estimations of the future COVID-2019 impact on the population, measures required from the public health system and effectiveness of different quarantine measures.

Recently, there is a significant number of models proposed to achieve this goal, varying from very simple to complicated ones, which include a lot of variables and parameters. However, it should be pointed out that detailed models, which include a lot of variable parameters require, respectively a detailed statistics for their validation. This is, for example, the case of Chinese outbreak, where the larger specification of affecting factors and decease statistics is available [5], [6], [7], [8] but for the rest of the world, the available data variety is quite limited. The last fact argues in favour of simple models, in majority regressions and compartmental ones. Among the last, one can mention variants of SEIR and SIR(D) models. e.g. [9], [10], [11], [12]. At the same time, even such models considered directly in their original differential form contain more parameters than directly available information for their determining that argues that the reduction of variable and parameters number may be preferable, see the discussion in [12]. On the other hand, the pure logistic growth fitting applied ad hoc simply as a probable regression for the expected sigmoid curve described the total number of infected persons leads to controversial and uncertain results, e.g. [13], [14], [15], that indicates an oversimplification when the fitting model does not have a clear linking to the epidemiological model background.

Thus, the main goal of this work is to fill a gap between two mentioned approaches trying to simplify the model as much as possible but not simpler and to check its reliability with respect to the actual available epidemiological data on COVID-19 outbreak and perspectives of its forecasting for the nearest months.

As a starting simplest model, the SIR model proposed by Kermack and McKendrick [16] is chosen. This approach, which is one of the most popular basic models of mathematical epidemiology [17] operates with the subdivision of the full population into three groups: Susceptible S, Infected I, and Recovered/Removed R, which interact according to the kinetic scheme where the rate coefficients are: k is the contact rate, or, in a general case the probability of being infected (this more abstract definition may be more reasonable when one considers not a localised well-mixed population but the dynamics of data averaged over a wide geographic area, see the respective discussion in [18]); τ is the characteristic duration of illness.

The system of ordinary differential equation, which corresponds to the kinetic scheme (1)–(2) reads as

Why may this simple system be appropriate for estimating the dynamics of the recent COVID-19 outbreak in different countries? The answer is based on the specificity of available data provided, e.g. by the World Health Organisation [19] and the European Centre for Disease Prevention and Control [20].

The freely available statistics is externally averaged since the daily data are given for each country as a whole, without more precise localisation within country regions and cities; detailed data on intra-country traffic, movements and contacts are absent. Therefore, the consideration can be based on the strong mean field level of approximation only if do not intend to introduce some hypothetical “hidden” parameters but intends to apply these available data only. In addition, the widespread restrictions of the last two months on international flights support the consideration of the country-based data.

On the other hand, it is not reflected in the available data, at what stage of illness each infected person enters the statistics as a registered infected one. The fact can be revealed from the latent phase during some testing of population groups suspicious to the risk to be infected up to the late phases of active illness when a person is taken to hospital under serious acute conditions. Whence, the division into latent (exposed) and active forms, i.e any kind of SEIR-models will be over-complexification with respect to the actual data uncertainty.

Thus, this model taken as a background will be a subject of study with respect to the further reducing to the form, which should provide a simple and easily realisable procedure as well as for testing of its relevance for modelling COVID-19 in different countries of the world.

Methods

Analytic treatment

The main fact, which allows the further consideration, is based on the proof given in [21] that the system of SIR Eqs. (3)– (3) is already one-dimensional due to the presence of two conservation laws and, moreover, can be sequentially reduced to the Verhulst (logistic) equation under the very natural assumptions. For the self-consistency of explanations, this line of reasoning is given below as well as its adjustment to the considered situation.

For simplicity, it is convenient to make this derivation for S, I, R considered as fractions of the full initially susceptible population taken as unity: . This gives the first conservation law reducing the number of equations from three to two.

Further, Eqs. (3) and (3) give with the solution written under the condition (this means that one practically neglect by the fraction of initially infected persons in the case of a large population). This gives the second conservation law that reduces the SIR system to one differential equation, which is convenient to write with respect to R:.

Under conditions of “weak outbreak” that means kτmax(R) <  < 1, or, if to consider the dimensional population where N is the total population, kτmax(R) <  < N (this, fortunately, fulfils for all countries, the cumulative number of infected is much less (e.g. for China, where the outbreak is almost finalized) and expected to be much less for the rest of world that the total country population), the exponential term in Eq. (3) can be expanded into the Taylor series that results in the classic Verhulst (logistic) equation where is the basic reproductive number.

Note that the transition to the dimensional quantities does not change its functional form and the initial growth rate () but simply rescales the saturation values from to .

Note also that the number of recovered/removed persons modelled by the Verhulst equation can be deduced from the number of infected persons I(t) (in the original dimensional form of data), due to Eq. (3) as the integral with varying upper limit

Note that Eq. (5) contains the characteristic disease duration τ, which is, in principle, known from medical practice but is not required from the pure data-driven approach to the problem. Actually, since τ is simple scaling constant, the replacement of R(t) by the cumulative number of infected persons does not change the form of the Verhulst equations written now as where r is the growth rate and K is the saturation value.

Its solution depends on a unique initial condition . However, the practical impossibility to determine the actual zero time of the outbreak, when the very first patient was infected, leads to the more convenient for data analysis purposes form of the solution of Eq. (7), where t defines the time moment when the cumulative number of infected persons reaches half of the saturated values. For the Verhulst equation, this moment of time also exactly corresponds to the maximum of the outbreak that is also practically important information from the prognostic point of view.

Numerical treatment

The data for analysis were downloaded from the official information page of the European Centre for Disease Prevention and Control [20] as a table in .csv format. They cover today’s data on the geographic distribution of COVID-19 cases worldwide starting from December 31, 2019 (in some cases, from the date of the first reporting to the ECDC database) up to the day of downloading, daily. The last day in the database used in the section “Results” is April 3, 2002.

The downloaded table was processed: for each country used for analysis identified by its geoID the column of daily infected cases was extracted and sorted in the descended order of dates (in the ECDC database, new data are added, conversely, to the top of the column). The first day of the records was identified and the time variable stated from this day as with daily step was formed. Such a two-columns matrix was saved as a file in .csv format for further parameter determination. The MATLAB script, which evaluated these procedures, is given in the Appendix.

The obtained data dale was processed using Loglet Lab 4 software (https://logletlab.com/) especially developed for high-quality fitting logistic and related models. A genetic algorithm (GA) based on Monte-Carlo Annealing was applied with the root mean square objective function; the uncertainty in parameter values was estimated using the bootstrap method. The Monte Carlo ensemble consisted of 104 iterations for both procedures, the confidence interval is equal to 95 %.

As one of the specific features applied in Loglet Lab for both fitting procedure and output representation, is the Fisher-Pry transform [22], which introduces the variable alteration (in the notation used in this work and in figures below) and transformation of the logistic function (8) transforms to the linear dependence

In addition to computational issue, this representation provides visual guidance for estimating time intervals, when the data already satisfies (if satisfies) the solution of the SIR model in the logistic representation that is one of the targets of this work.

Results

The first test evaluated is naturally related to China, where the active outbreak is finished, and the situation with the cumulative number of infected persons reached the saturation, see Fig. 1 (left panel). For this quantity, the fitting with the solution (8) not only visually gives reproducing the actual dynamics but also has a relatively high value for the coefficient of determination . However, ΣI is an integral value, which smooths instant dynamics for the daily number of infected persons.

Fig. 1

The comparison of COVID-19 outbreak in China (markers) and its fitting with the SIR model (solid line) reduced to the logistic equation with respect to ΣI (left subpanel) and other characteristics: daily I (rate, middle panel), and the Fisher-Pry representation of the logistic growth, where .

This takes place actually, see Fig. 1(middle panel), where the real data markers correspond to the solid line curve representing the rate of change of the logistic curve more qualitatively that quantitatively.

In more details, one can note that they go close to each other during the first stage of the outbreak but start deviate drastically not reaching the maximum, which is completely “cut-off” the peak and lead to the fast decay of the number of infected persons. This effect is definitely connected with the quarantine measures introduced by the Chinese government that is supported by several studies building more complex models aligned with the timeline of these measures, e.g. [6], [7].

At the same time, the last stage of the decaying outbreak again follows the bell-shaped curve given by Eq. (9) that is more-or-less natural from the point of view of the theory of dynamical systems since the solution goes fast to the stationary point via the universal exponential function.

As well, the total duration of the outbreak is also well reproduced by the simple logistic model. It is confirmed by Fig. 1(right panel), where the Fisher-Pry plot is shown. Although the real data markers visibly fluctuate around the fitting straight line, they are most close to the latter in the first growing stage of the outbreak, i.e. the prognosis of the total duration made before the expected maximum gave the correct total duration. Thus, the main effect of quarantine measures, which changed the functional form of the plot was mainly in a drastic reduction of the number of affected people during these predicted time interval that is definitely the most important form the point of view of public health.

However, the current situation in several other countries argues in favour of better applicability of the proposed simplified model. Italy can be considered as the most demonstrable example since it is a country the most heavily affected among European ones now and the time course of the outbreak seems to leave the maximum behind. This provides an opportunity to explore a predictive capacity of the simplified model taking into consideration different time intervals from the outbreak’s beginning as shown in Fig. 2 . The respective of the parameters and their uncertainty found from the data fitting with different lengths of time series are given in Table 1 .

Fig. 2

Series of plots, which demonstrate logistic curve fittings (solid parts of curves) and predictions(dashed parts) using the data (markers) registered for different time intervals from the outbreak’s beginning (80, 85, 90, and 95 days since 31/12/2019, day 0). The yellow filled areas denote uncertainty ranges obtained via the Monte Carlo simulations. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Table 1

Parameters of the data fitting for Italy corresponding to Fig. 2. Here Fitted means the direct fitting of the data, Min, Max, and Median are the boundary values and the median for the distribution obtained via the Monte Carlo simulations (note that the order of magnitudes for r values is reverse to the order for K and t but they are kept in this order to keep the set of parameters for each curve in the same column. Dimensionalities for K, r, and t are Persons, day and day, respectively. The notation in parenthesis denotes the number of days since 31/12/2019, for which the data were used for fitting.

Parameter	Fitted	Min	Max	Median
K (80 days)	84784	70917	124125	92812
r (80 days)	0.209	0.219	0.190	0.203
tm (80 days)	81	79	84	81
K (85 days)	123148	113506	136838	123024
r (85 days)	0.194	0.201	0.186	0.194
tm (85 days)	83	83	84	83
K (90 days)	121848	118000	135206	122401
r (90 days)	0.193	0.197	0.179	0.191
tm (90 days)	84	83	85	84
K (95 days)	132455	130766	152018	133371
r (95 days)	0.183	0.187	0.161	0.182
tm (95 days)	85	85	87	85

From the set of plots in Fig. 2 one can see that in the case of epidemic outbreak in Italy, in contrast to the case of China, the correspondence between the solutions of the SIR model reduced to the logistic equation and the real recordings is quantitatively accurate. It is seen not only from the plots of the cumulative number ΣI ( for the first plot based on the shortest sample, and for the rest three) but also in the subpanels illustrating the daily number of infected persons, i.e. the rate of change of the logistic function, where markers are slightly and symmetrically scattered around the bell-shaped functions (9) as well as in the Fisher-Pry plots, where the data accurately follow the straight line after some short setting interval affected by small numbers of the daily data counts.

Thus, it is possible to conclude that the SIR-model in its logistic growth representation allows the accurate fitting the epidemic process in the case of Italy from the accelerated growth stage to the stage of the epidemic’s decay (decelerated growth).

However, another important question arises: what is the capacity of this model for predicting further epidemic course. The sequence of fittings limited by different final dates separated by ten-days intervals gives information for such analysis. First of all, it should be pointed out that the growth of sample leads to the diminishing value of the coefficient r and the subsequent growth of the coefficient K (the maximal estimated number of infected person that is on of the key demanded characteristics expected from the prognosis). Such behaviour is not only typical for the logistic regression [23] but also has an explanation based on the model’s constructions: comparing Eqs. (4) and (7) one can see that both r and K depend on one parameter ℜ0 reciprocally. Whence, the decay of r from 0.209 down to 0.183 is not connected with the actual reducing the averaged contact rate (this case would result in the changing shape of the rate’s bell curve as in the case of China that is not observed in the case of Italy) but with the sequential better and better adjustment of the logistic curve to the longer and longer time series of data. This conclusion is also supported by the fact that all Fisher–Pry plots are accurately linear.

At the same time, an application of the Monte Carlo method for estimating parameters uncertainty provides some opportunity for forecasting and estimation its reliability even when the outbreak still does not reach its peak. This method not only gives a range of predicted values, which is wider for less warrant early forecast and squeezes to the purely fitted value when the processed sample gets longer, but also propose some criteria for estimating relevance of the extrapolation. In particular, there is visible trend in the displacement of the predicting extrapolation of the simple fitting to the lower boundary of the prediction uncertainty defined by the Monte Carlo method. An explanation of this trend is based on the better relevance of the fitting for the data, which are in the regime of decelerating growth, i.e. the rate maximum is behind. Thus, the uncertainty in the data, which will be added, lead to the possible growth of the cumulative variable considered. On the contrary, before the rate maximum both uncertain results are possible, over-fitting as well as under-fitting. From the practical point of view, for predictions on earlier stages, it is better to take into account the Monte Carlo-based upper boundary, which as it seen in Fig. 2 and Table 1 gives reasonable estimations even with such a simple model is used.

To determine, either the case of Italy is an exception or the reduced SIR model is already wide applicable, such tests were evaluated for a number of other countries and argued in the positive answer for the majority of them. Some typical examples for different stages of outbreak are shown in Fig. 3 , where the Fisher-Pry representations on the date of April 3rd of 2020 are given.

Fig. 3

The Fisher-Pry plots for different countries based on the data available to April 3, 2020 (markers). The solid and dashed line show their fitting and forecast, respectively; the yellow area covers the uncertainty range. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

One can see that five out of six given examples follow the linearity of the Fisher–Pry plot; moreover, this behaviour practically does not depend on the current stage of the epidemic outbreak and is observed for Germany and Spain, where the raw data and fitting indicate the decelerating phase, for France, which seems to be close to the peak as well as for Russia where the data argue the first, accelerating, phase. Except for France, where the three other countries indicate similar r ≈ 0.215 ± 0.015 by the line’s slope with uncertainties diminishes as the data approaching the latter phase. As for the USA, the Fisher-Pry plot is also linear but with the drastically larger . However, it is clearly seen that the SIR is not universal for the COVID-19 pandemic. In addition to the partial deviation for China visible in Fig. 1, Fig. 3 indicates that South Korea is a complete exception: the plot in Fisher–Pry’s variables is not linear almost everywhere; the only short linear interval is located around the middle point (50%), where the maximum of I occurs, i.e. is simply connected with the linearity of ΣI around this point, therefore it is a pure mathematical artefact.

Discussion and conclusions

The main conclusion is that the evaluated analysis shows that the SIR being sequentially reduced to the Verhults equation can provide an accurate description of the COVID-19 epidemic on the level of a country-wide averaged statistics. Moreover, the time course of outbreaks follows such functional dependence practically as a rule. The clearly visible exceptions are exhibited by the countries, where very strong and wide-ranged measures were undertaken. First of all it is China, and the change of epidemic dynamics is visible in Fig. 1(middle panel) that allows to assign and localise in time an effect of this measures be a deviation from the “natural logistic course” typical for the majority of countries. The second demonstrable example is South Korea, where completely different but also effective kind of measures was undertaken [24]. This is also reflected in the respective subpanel of Fig. 3, where the Fisher-Pry representation of the data differs from both China and other countries represented there. Thus, it is possible to conclude that the proposed SIR-Verhults representation can play the role of a criterion for quantification such public health measures.

For other countries, for which illustrations are given in Figs. 2 and 3 as well as many others the strictly logistic growth of the cumulative number of infected persons is demonstrated as typical. Note also a recent study [15] the logistic regression for a wide range of countries and their territories was evaluated. However, the solution of the Verhulst equation was applied just as a kind of sigmoidal function without its linking to any epidemiological model, and, as a result, reliability of the obtained parameters were discussed on the base of pure computational issues only.

On the contrary, in this work, the parameters of the logistic function are proven as connected with the parameters of the well-established Kermack–McKendrick model that reveal their meaning in the context of the epidemic process. In particular, the proven interrelation between r and K parameters indicates that the diminishing of r during the outbreak can originate not from reducing contact rate but principally from the better regression accuracy when the time series get longer and take into account the normal effect of deceleration characterising the full interval of the outbreak. Simultaneously, this provides an opportunity to estimate the quality of the logistic curve-based forecast and its uncertainty more warrantly.

Form this point of view, the Fisher-Pry plot is significant since it indicates the logistic character of the cumulative number of infected persons when this plot is linear even if its slope determined for time series of different length is different. This is a reason why only such plots are given for different countries in Fig. 3. In addition, the length of data series grows daily; therefore more and more accurate estimations will be available but their possibility is assured namely by this form of data and model representation.

Declaration of Competing Interest

The author declare that he has no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.