Power-law multi-wave model for COVID-19 propagation in countries with nonuniform population density

Abstract

Background: The purpose of our study is to develop a quite precise mathematical model which describes epidemics spread in a country with non-uniform population density. This model gives explanation of quite long duration of the peak of a respiratory infection such as the coronavirus disease 2019 (COVID-19).

Methods: The theory of kinetic equations and fractal analysis are used in our mathematical model. According to our model, COVID-19 spreading takes the form of several spatio-temporal waves developing almost independently and simultaneously in areas with different population density. The intensity of each wave is described by a power-law dependence. The parameters of the dependence are determined by real statistical data at the initial stage of the disease spread.

Results: The results of the model simulation were verified using statistical data for the Republic of Belarus. Based on the developed model, a forecast calculation was made at the end of May, 2020. It was shown that the epidemiological situation in the Republic of Belarus is well described by three waves, which spread respectively in large cities with the highest population density (the first wave), in medium-sized cities with a population of 50−200 thousands people (the second wave), in small towns and rural areas (the third wave). It was shown that a new wave inside a subpopulation with a lower density was born 20−25 days after the appearance of the previous wave. Comparison with actual data several months later showed that the accuracy of forecasting the total number of cases for a period of 3 months for total population in the proposed approach was approximately 3%.

Conclusions: The high accuracy mathematical model is proposed. It describes the development of a respiratory epidemic in a country non-uniform population density without quarantine. The model is useful for predicting the development of possible epidemics in the future. An accurate forecast allows to correctly allocating available resources to effectively withstand the epidemic.

Graphical abstract

Keywords

COVID-19 / forecast model / simultaneous waves / population density

1 INTRODUCTION

Predicting the spread of the coronavirus disease 2019 (COVID-19) is one of the most challenging problems for the relevant scientific communities of medics, mathematicians and physicists [1,2]. An accurate forecast allows to correctly allocate available resources to effectively withstand the epidemic [3,4]. It is impossible to predict the spreading of COVID-19 with absolute accuracy. However, forecasting can be improved by adjusting the mathematical models used and comparing the forecasting data with the actual situation.

A characteristic feature of the development of the epidemic situation in some of countries was a long plateau in the number of new cases of the COVID-19 or a slow decrease in daily cases in April−May 2020 [5]. This was observed in the United States (April and May), United Kingdom and Canada (from the beginning of April to the middle of May), Russia (in May), Belarus (from the last decade of April to the last decade of May), Sweden (April and May), Indonesia (April and May), Poland (April and May), Ukraine (April and May). To note that this situation with long plateau is characteristic for the total epidemic indicators worldwide, when the total number of sick people is at the level of 80−100 thousand new cases per day starting from the first decade of April to the last decade of May.

For fast diseases, with an incubation period of several days [6], the presence of a long plateau is not typical situation. More typical situation is the epidemic development according to the scenario “growth−peak phase−decline”, when the duration of the peak phase is several times shorter than the growth phase [7]. For the cases under discussion, the plateau duration is comparable and even exceeds the duration of the rapid growth phase (approx. 30 days and more). This feature raises questions about our understanding of this phenomenon.

To the middle of spring 2020, there were reports about the expectation of the second wave of the COVID-19 spread [7]. The development of the second epidemic wave in many countries in the fall of 2020 is no longer in doubt. Note that this refers to the second wave of the disease in the same population in some time. But we consider the phenomenon of the simultaneous propagation of several waves in different subpopulations of the same population. This phenomenon can explain the current situations with long plateau. Below is presented the mathematical model that describes this phenomenon.

2 MATHEMATICAL MODEL

As a result of the analysis of the possible reasons of the long plateau phenomenon, we came to the conclusion that it can be associated with regional characteristics of the nonuniform distribution of the population density and the formation of the additional waves of the virus decease spreading in corresponding subpopulations. To note that in this work we do not consider well known effect of the presence of a large number of asymptomatic cases of a mild course of the COVID-19 [7]. In fact, such cases are not amenable to counting and for analysis in real epidemic situation. We will only use information about reported cases, which are usually associated with more severe manifestations of the disease. Using the open statistics for the Republic of Belarus [5], we want to show that a long plateau in this case can be described by the simultaneous propagation of three independent waves of COVID-19 in the country.

The essence of this phenomenon, in our opinion, is as follows. At the beginning of the epidemic, all countries closed their borders. Therefore, each of the countries, including the Republic of Belarus, can be considered as a closed and isolated human community. However, within many countries the population density is nonuniform one. It can be considered as quite reliable fact that a population density is one of the key factors influencing the spread of the respiratory virus [8,9]. In large cities with a high density, the spread of the virus is most intense.

From the point of view of population density, the Republic of Belarus can be divided into 3 agglomerations: 1) the Minsk city with the population density of 5800 ppl km⁻²; 2) regional centers, where the population density is from 2400 ppl km⁻² (Brest), 2600 ppl km⁻² (Grodno) to 3800 ppl km⁻² (Gomel), and large cities of regional subordination (Polotsk, Orsha, Bobruisk, Pinsk, etc.), with a population density of 2100 ppl km⁻² (Polotsk) to 2950 ppl km⁻² (Orsha); 3) medium and small district centers and rural settlements with a population density of tens to hundreds of people per square kilometer.

It is worthy to note that the average population density for the Republic of Belarus is 46 ppl km⁻². The spread of the virus in each of these agglomerations should occur at its own speed, the onset of the disease in each of the agglomerations can also be shifted in time [8]. We assume that the spread of the disease begins in the largest city with the highest population density and then, with a certain delay, the disease comes in an agglomeration with a lower density. In the first approximation, it can be considered that the disease in each of the agglomerations proceeds independently and mainly propagates within a subpopulation. These considerations underlie the analysis below.

Previously, we used the described mathematical approach [10] to analyze and forecast the situation with the dynamics of the coronavirus spreading in the Republic of Belarus [11–13]. The essence of the analysis is to graph the dependence of new cases of the disease on cases in the previous day in double logarithmic coordinates (Fig.1) [11]. The authors of paper [10] found a general pattern for a number of countries, namely, the power-law dependence of the number of new cases on a given day n+1 on the number of cases on the previous one:

(1)C(n+ 1)=\alpha C(n{)}^{\beta }.

Here n is the days counted from the beginning of the spread of the epidemic disease. The observation of similar power-law dependences is reported also in other, earlier works [14]. The reason for the good description of the epidemic situation by Eq. (1) is not fully understood. Perhaps this is due to the heterogeneity of the human community, manifested in both the heterogeneous population density [15,16] and age and activity heterogeneities [17]. This requires additional research.

Since this relationship underlies our mathematical consideration, we called the proposed model as phenomenological one. The exponent \beta in Eq. (1) and the proportionality coefficient \alpha should remain constant over time if the intensity of social contacts and the scenarios of social behavior in the population did not change. The presence of such dependence makes it possible to determine its parameters at the initial stage of the epidemic and to forecast a change in the disease rate over quite long time periods.

It can be shown [10] that the Eq. (1) leads to a geometric progression, which in turn gives an analytical expression for the number of diseased people C(n) over time, i.e. from the number of days n from the start of the count:

(2)C(n)={\alpha }^{\frac{1-{\beta }^n}{1-\beta }}C(0{)}^{{\beta }^n}.

To build a forecast curve, one also needs the initial value of the number of involved people C(0). It is usually accepted as the starting day to use the date when the total number of cases in the country has reached a certain level. Since the disease statistics always show significant fluctuations at the initial stage, Eq. (1) is applicable when the number of cases is about 100 people or more. Initially the parameter C(0) in our work was determined from the condition of the minimum dispersion in approximation of the actual data in the period from March 30 to April 24:

(3)C(0)\approx 100.

As noted above, the day when the total number of cases exceeded 100 is selected as the starting point for our simulation. For the Republic of Belarus, this is March 30, 2020.

The considered mathematical approach is in agreement with the logistic equation of the spreading of virus diseases [18–20]. The Eq. (1) for the dynamics of the viral disease spread is closely related to the classical approach based on the Verhlust logistic equation. In similar notation, this equation can be written as [21]

(4)\frac{dC(t)}{dt}=kbC( t)\left\{1-\frac{C(t)}{N}\right\}.

Here k is the probability of transmission of the infection from the patient to a healthy person, averaged over the entire population, b is the average number of contacts per day for the average person, N is the population size (population of the country), t is temporal variable. Note that the well-known viral disease reproductive number R_0 [22] is written in these notation as R_0=kb. The solution of Eq. (4) with the initial condition C{|}_{t=0}= C(0) is well known and has the form

(5)C(t)=\frac{NC(0)\exp (R_{0}t)}{N+C(0)[\exp (R_{0}t)-1)]}.

Then, for short times, t\ll t_{\ast }=1/R_0, the number of cases should grow exponentially

(6)C(t)\approx C(0) \exp (R_{0}t).

Comparing Eqs. (2) and (6) and taking into account that the small parameter is 1-\beta \approx 0.05, we can establish an approximate relation between the parameters

(7)R_0=k b\approx \ln \alpha +(\beta -1)\ln C(0) ,t,n\ll t_{\ast }.

According to the Eq. (7), an increasing in values of the parameters \alpha and \beta indicate a higher reproduction rate in the corresponding population.

Additional data accumulated to mid-May showed a certain deviation from the forecast curve based on the single wave model [11,12] (Eqs. (1), (2); Fig.2). Initially, we assumed that the observed morbidity situation can be described more precisely if its propagation is allowed in the form of two waves. These waves independently but simultaneously arise, develop and propagate in the two largest population agglomerations described above. The addition of the second wave to the model well improved the compliance between the calculation results and the actual data in the period from 40 to 50 days from the start of the count. But then again, a certain deviation occurred in the daily cases of the COVID-19 disease (Fig.3). This suggested that, nevertheless, we are dealing with three simultaneous and independent waves in the three subpopulations described above [13].

The described nature of the development of three waves of COVID-19 can be confirmed by their initial time shift relative to each other. To obtain the numerical parameters of the model, the following algorithm was proposed.

The characteristics of the first, most intense, wave can be described by parameters obtained from data at the initial stage of the development of the epidemic (March 6 − April 30). These were done in the preprint [11]. New data and their deviation from the first forecast curve for one propagation wave can be used to identify the parameters of the second wave. For this, data from 33 to 50 days were used [12] (Fig.3). Finally, additional data from 51 to 57 days were used to identify the parameters of the third wave [13].

Thus, we put forward the hypothesis that the epidemic situation with COVID-19 in the Republic of Belarus can be described by three independent waves:

(8)C(n)=C_1(n)+C_2(n)+C_3(n).

For the first wave, according to the above assumptions, the following relations are true:

(9)C_1(n+1)={\alpha }_1{C}_1(n{)}^{{\beta }_1}.

(10)C_1(n)= {\alpha }_1^{\frac{1-{\beta }_1^n}{1-{\beta }_1}}{C}_1( 0{)}^{{\beta }_1^n}.

The second wave is described by the similar relationships with new constants:

(11)C_2(n_2+ 1)={\alpha }_2{C}_2(n_2{)}^{{\beta }_2},

(12)C_2(n_2)={\alpha }_2^{\frac{1-{\beta }_2^{n_2}}{1-{\beta }_2}}{C}_2(0{)}^{{\beta }_2^{n_2}},

(13)n_2=n -\mathrm{\Delta }{n}_2.

Finally, the third wave is described by similar relations:

(14)C_3(n_3+ 1)={\alpha }_3{C}_3(n_3{)}^{{\beta }_3},

(15)C_3(n_3)={\alpha }_3^{\frac{1-{\beta }_3^{n_3}}{1-{\beta }_3}}{C}_3(0{)}^{{\beta }_3^{n_3}},

(16)n_3=n -\mathrm{\Delta }{n}_3.

Here \mathrm{\Delta }{n}_2, \mathrm{\Delta }{n}_3 is the delay time of the second and third waves relative to the first (in days).

3 RESULTS AND DISCUSSION

Thus, according to the available data and proposed model (8)−(16), it is necessary to determine 11 parameters for the first, second and third waves: {\alpha }_1, {\beta }_1, C_1(0), {\alpha }_2, {\beta }_2, C_2(0), \mathrm{\Delta } n_2 and {\alpha }_3, {\beta }_3, C_3(0), \mathrm{\Delta } n_3. As the objective function, the standard deviation of the calculated data and actual ones on the total number of cases for the entire period was considered.

The search for optimal parameters was carried out using the gradient method of the steepest descent [23]. We found the minimum of the total deviation of the actual data from the calculation data. The results of the search are:

(17){\alpha }_1=1.72(4), {\beta }_1=0.94(9), C_1(0)= 100,

(18){\alpha }_2=1.72(7), {\beta }_2=0.94(2), C_2(0)= 24,\mathrm{\Delta }{n}_2=25,

(19){\alpha }_3=1.72(7), {\beta }_3=0.94(2), C_3(0)= 9,\mathrm{\Delta }{n}_3=39.

It is important to note that the second and third waves are described by identical parameters \alpha and \beta, which are similar in magnitude to the parameters of the first wave. The difference is only in the third decimal place. All parameters coincide within the error of the data processing method used. The second and third waves differ, as we expected, with a lower initial value of the incident cases C_2(0)=24, C_3(0)=9, and are characterized by a shift of approximately 25 days and 39 days relative to the first wave. In a first approximation, it can be assumed that the internal patterns of the distribution of the three epidemic waves are approximately the same (parameters \alpha and \beta). Waves differ only in their intensity (amplitude) and onset time.

Close values of the parameters \alpha and \beta for all three waves were obtained as a result of mathematical processing of the primary data. It can be noted that in [10] close values of these coefficients were obtained even for three countries. And these coefficients are close to our values. We do not have an exact explanation on this. It can be assumed that this is a manifestation of the fact that in different places of the same country the same cultural traditions and scenarios of social behavior are adopted (shaking hands; getting close or far when talking; gathering on weekends in companies, etc.).

The results of the predictive calculation in the three-wave approximation [13] by the Eq. (8)−(19) are presented in Fig.4. The error in our calculations is preliminary estimated at 10%−15% at time intervals of about 1 months. The data shown in Fig.2, Fig.3 indicate the correctness of such estimates. In fact the forecast accuracy is about 3% for the total number of cases and no more than 15% for daily new cases within 3 months from the date of the forecast. The performed forecast calculations are valid as long as the conditions for the development of the morbidity process correspond to the conditions under which the initial statistics for the forecast were collected.

For a better understanding of the accuracy of various approximations, we presented the results of consideration in the approximation of one, two, and three waves (Fig.2-Fig.4, respectively). For clarity, the contribution of each wave is highlighted in color. As can be seen, the most accurate is the approximation of the three-wave model (Fig.4). From the results obtained above, we can conclude that the puzzle of the long plateau and general behavior in COVID-19 dynamics is explained by our phenomenological model. But the construction of such a model requires knowledge of the regional characteristics of the country and details of its population density.

Let’s briefly discuss the issue of the accuracy of forecast models for the spread of epidemics. A special article was published on the situation with COVID-19 with discussion the problem of low forecast accuracy [24]. In many cases, the error in the forecasts is measured not in percent, but hundreds percent and even orders of magnitude. The reason of such discrepancy is poor data input, wrong modeling assumptions, high sensitivity of estimates, consideration of only one or a few dimensions of the problem at hand, selective reporting etc. Indicative in terms of accuracy is paper [25]. In this paper similar forecast have been made about the propagation of the coronavirus epidemic at the initial stage of spread in various states of India. The model used in this work contains 6 free parameters that must be found from the analysis of primary data. In this case, a period of about 40 days is considered. From the data of this work, it can be seen that the forecasting error reaches 20%−30% in some cases. In another, later work written by the same authors, an updated short-term forecast for a period of 20 days is constructed ([26], Fig. 7). The error of such a forecast reaches 100% or more for states with a small number of daily cases (several dozen per day). But for a state with a large number of new daily cases (about 500 per day), the short-term forecast error is approx. 10%. From this point of view, the model proposed in this work can be considered as sufficiently accurate.

4 CONCLUSIONS

The phenomenological model is proposed that describes the dynamics of the spread of COVID-19 in Belarus, a country with population about 10 million people. We used the approximation of the propagation of three independent waves in subpopulations with different population densities. The dynamics of each of the waves is described by a universal power-law dependence, first discovered in Ref. [10]. The forecast calculations performed using the proposed model on the base of statistical data for the Republic of Belarus showed a high forecast accuracy. The model is based on limited statistical information on daily morbidity from March 6 to May 25 (total 81 points, 57 points were used for analysis). Nevertheless, the forecast accuracy is about 3% for the total number of cases and no more than 15% for daily new cases within 3 months from the date of the forecast.

We hope that the proposed approaches will be useful for analysis of the pandemic situation in other countries. It is possible that for large countries where the distance between large cities exceeds 500−600 km (the average distance available for a one-day trip by car), the propagation dynamics of COVID-19 also will be described by the superposition of a several simultaneous spatio-temporal waves. The proposed approach can be applied also for analysis of a spread of seasonal respiratory viral infections, such as influenza.

References

[1] AnastassopoulouC.,RussoL.,TsakrisA.. (2020) Data-based analysis, modelling and forecasting of the COVID-19 outbreak. PLoS One, 15, e0230405

[2] PetropoulosF.. (2020) Forecasting the novel coronavirus COVID-19. PLoS One, 15, e0231236

[3] Grinchuk,P. S., Fisenko,K. I., Fisenko,S. P. Danilova-Tretiak,S. (2021). Isothermal evaporation rate of deposited liquid aerosols and the SARS-CoV-2 coronavirus survival. Aerosol Air Qual. Res., 21 : 200428

[4] SiettosC. I.. (2013) Mathematical modeling of infectious disease dynamics. Virulence, 4, 295–306

[5] COVID-19 Dashboard by the Center for Systems Science and Engineering at Johns Hopkins University.

[6] LauerS. A.,Grantz K. H.,BiQ.,JonesF. K.,ZhengQ., MeredithH. R.,AzmanA. S.,ReichN. G.. (2020) The incubation period of coronavirus disease 2019 (COVID-19) from publicly reported confirmed cases: estimation and application. Ann. Intern. Med., 172, 577–582

[7] GlezenW. P.,DeckerM.,JosephS. W. MercreadyR. G.. (1987) Acute respiratory disease associated with influenza epidemics in Houston, 1981‒1983. J. Infect. Dis., 155, 1119–1126

[8] PedrosaR. H.. (2020) The dynamics of COVID-19: weather, demographics and infection timeline. medRxiv,

[9] Rocklöv,J and Sjödin, H. (2020) High population densities catalyse the spread of COVID-19. J. Travel Med., 27, taaa038

[10] FanelliD.. (2020) Analysis and forecast of COVID-19 spreading in China, Italy and France. Chaos Solitons Fractals, 134, 109761

[11] GrinchukP.. (2020) The forecast of the dynamics of COVID-19 in the Republic of Belarus. Adjusted data, Preprints. ru 1–11 (in Russian)

[12] GrinchukP. S.FisenkoS.. (2020) Two waves in the COVID-19 dynamics in the Republic of Belarus, Preprints. ru 1–14 (in Russian)

[13] GrinchukP. S.FisenkoS.. (2020) Several simultaneous waves of COVID-19. Explanation of long-time plateau, Preprints. ru 1–15 (in Russian)

[14] StroudP. D.,SydoriakS. J.,RieseJ. M., SmithJ. P.,Mniszewski S. M.RomeroP.. (2006) Semi-empirical power-law scaling of new infection rate to model epidemic dynamics with inhomogeneous mixing. Math. Biosci., 203, 301–318

[15] NeipelJ.,Bauermann J.,BoS.,HarmonT.. (2020) Power-law population heterogeneity governs epidemic waves. PLoS One, 15, e0239678

[16] CampiG.,Mazziotti M. V.,VallettaA.,RavagnanG.,MarcelliA., PeraliA.. (2021) Metastable states in plateaus and multi-wave epidemic dynamics of COVID-19 spreading in Italy. Sci. Rep., 11, 12412

[17] BrittonT.,Ball F.. (2020) A mathematical model reveals the influence of population heterogeneity on herd immunity to SARS-CoV-2. Science, 369, 846–849

[18] BlumbergA.. (1968) Logistric growth rate functions. J. Theor. Biol., 21, 42–44

[19] TsoularisA.. (2002) Analysis of logistic growth models. Math. Biosci., 179, 21–55

[20] VerhulstP.. (1838) Notice sur la loi que la population poursuit dans son accroissement. Corresp. Math. Phys., 10, 113–121

[21] NahushevA.. (1995) Equations of Mathematical Biology (Higher School, Moscow)

[22] LiuY.,Gayle A. A.,Wilder-SmithA.. (2020) The reproductive number of COVID-19 is higher compared to SARS coronavirus. J. Travel Med., 27, taaa021

[23] KornG. A.Korn T.. (2000) Mathematical Handbook for Scientists and Engineers: Definitions, Theorems, and Formulas for Reference and Review (Dover Civil and Mechanical Engineering). Dover Publications

[24] IoannidisJ. P. A.,CrippsS.TannerM.. (2022) Forecasting for COVID-19 has failed. Int. J. Forecast., 38, 423–438

[25] Sarkar,K., Khajanchi,S., Nieto,J. (2020). Modeling and forecasting the COVID-19 pandemic in India. Chaos Soliton. Fract., 139 : 110049

[26] Khajanchi,S., (2020). Forecasting the daily and cumulative number of cases for the COVID-19 pandemic in India. Chaos 30, 071101

RIGHTS & PERMISSIONS

The Author (s). Published by Higher Education Press.