Modeling the multifractal dynamics of COVID-19 pandemic.

To describe the COVID-19 pandemic, we propose to use a mathematical model of multifractal dynamics, which is alternative to other models and free of their shortcomings. It is based on the fractal properties of pandemics only and allows describing their time behavior using no hypotheses and assumptions about the structure of the disease process. The model is applied to describe the dynamics of the COVID-19 pandemic from day 1 to day 699 from the beginning of the pandemic. The calculated parameters of the model accurately determine the parameters of the trend and the large jump in daily diseases in this time interval. Within the framework of this model and finite-difference parametric nonlinear equations of the reduced SIR (Susceptible-Infected-Removed) model, the fractal dimensions of various segments of daily incidence in the world and variations in the main reproduction number of COVID-19 were calculated based on the data of COVID-19 world statistics.

Introduction

At present, one of the most demanded issues of mathematical modeling is developing adequate mathematical models of pandemics. The relevance of constructing and applying such models has been greatly enhanced by the COVID-19 pandemic.

The quantitative description of pandemics is based on the function v(t), which determines the time dependence of the number of daily diseases in the world associated with specific pandemics. This function is determined using relevant statistical data. In many papers, the function v(t) is determined by solving differential and integro-differential equations.

The function v(t) constructed basing on the world statistical data for COVID-19 pandemic determines a complex multifractal curve. As shown by calculations of the fractal dimensions D for various segments of v(t) during the time interval of interest for us, the values of D lie within an interval from 1.0708 to 1.4118. From this fact it follows that the curve v(t) is not differentiable. Therefore, the applicability of differential and integro-differential equations to the pandemics description becomes problematic.

The description of v(t) evolution is greatly complicated by the presence of jumps. For example, on December 10, 2020 a substantial jump of COVID-19 daily incidence in the world had been observed. On this day, the daily incidence increased by 2.4 times as compared to the preceding time. Then this number fell to preceding values during a day.

Within the approaches based on differential or integro-differential equations, the description of such jumps is not possible. Alternative descriptions of COVID-19 dynamics using fractional derivatives have been developed, e.g., in Refs. [1], [2], [3], [4], [5], [6], [7], [8], [9], [10], [11].

Based on the original Mandelbrot “Multifractal Walk down Wall Street” model [12], a model of multifractal dynamics (MFD) was developed and tested [13], [14], [15]. In this model, the fractal properties of the curve v(t) allow describing the pandemic dynamics using no assumptions and hypotheses about the structure of the disease process. In this case, the multifractality means that the entire temporal interval of observing the process is divided into subintervals, so that the fractal dimension has a certain value for each of these intervals. We believe that the extension of the multifractal approach to the COVID-19 pandemic dynamics will allow a fresh look at this topical problem and attract the attention of researchers in this field.

We start the present work basing on our experience in the field [13], [14], [15]. First, we calculate the fractal dimensions of the phase space of the COVID-19 pandemic and various segments of the daily incidence in the world according to the global COVID-19 statistics [16], [17]. Second, we estimate the variations in the basic reproduction number of COVID-19 [18] using the equations of the discrete delay model [20] or the reduced Susceptible-Infected-Removed (SIR) model [21], [22] to clarify short-term forecasts and identification of new trends in the dynamics of the COVID-19 pandemic [24], [28], developed within both the traditional and fundamentally new approaches and software tools [19], [22], [23], [25], [26], [27].

The paper is organised as follows. Section 2 gives the main definitions of multifractal dynamics with piecewise linear trends and the fractal dimension of the phase space of the COVID-19 pandemic. Section 3 simulates the multifractal dynamics of the COVID-19 pandemic. Section 4 is devoted to modeling the multifractal dynamics of the main reproduction number of COVID-19. Section 5 summarizes the results and perspectives. As initial data, we used the data of world statistics of COVID-19 [16] for the time interval from January 22, 2020 to July 6, 2021, indicated below as [1a] and [1b].

Multifractal dynamics with piecevise linear trends

The multifractal dynamics (MFD) model proposed to describe the dynamics of COVID-19 pandemic is based on fractal properties of the daily increment (rate) of the infected number curve, v(t) = ΔN(t) = N(t + 1) − N(t), where N(t) is the number of infected people at the moment of time t, the time t is measured in days. The total population is N max (for the world N max = 7.9 ⋅ 109). MFD allows describing the COVID-19 pandemic dynamics using the data of COVID-19 World Statistics without assumptions and hypotheses about the structure of the disease process, as shown in Fig. 1, Fig. 2 .

Fig. 1

The total number of coronavirus cases N(t) from the statistical data [1a] [16] or [1b] [17](left) and its approximation by the lower curve calculated from Eq. (3)(right).

Fig. 2

Daily new cases v(t) in units of N = 106 from the statistical data [1a] [16] (left) and [1b] [17] (right).

The dynamics of a multifractal process in the interval T : t ∈ [t ,  t ] can be divided into two components using the notion of a piecewise linear trend:where is the average rate, is the linear trend of a process that smoothly varies with timeand are fast oscillations relative to the trend.

The rate X of the linear trend has the meaning of a piecewise continuous acceleration (deceleration) of rectilinear uniformly variable motion over the plot of the functionwhich approximates N(t) in the interval T : t ∈ [t ,  t ].

As an example, we consider the statistical data [1a] and [1b] [16], where we divide the time interval 1 ≤ t ≤ 491 since January 22, 2020 till May 26, 2021 into ten intervals (with attached eleven intervals since May 26, 2021 till June 19, 2021 shown in Table 1 ) and construct ten piecewise linear trends , i = 1, 2, …, 10 on them. Fig. 1 illustrates the visual quality of approximations of the number of infected people N(t)/N by Eq. (3) and the daily increment of the infected number v(t)/N shown in Fig. 2 by means of piecewise linear trends of Eq. (2) displayed in Fig. 3, Fig. 4 .

Table 1

Dividing the interval from 22 January 2020 till 29 September 2021 into sub-intervals T1, …, T10∗, T11∗.

T1	[ 1, 50(50)]	[January 22, 2020  − March 11, 2020],
T2	[ 50, 75(25)]	[March 11, 2020  − April 5, 2020],
T3	[ 75,125(50)]	[April 5, 2020  − May 25, 2020],
T4	[125,165(40)]	[May 25, 2020  − July 4, 2020],
T5	[165,255(90)]	[July 4, 2020  − October 2, 2020],
T6	[255,285(30)]	[October 2, 2020  − November 1, 2020],
T7	[285,355(70)]	[November 1, 2021  − January 10, 2021],
T8	[355,390(35)]	[January 10, 2021  − February 14, 2021],
T9	[390,460(70)]	[February 14, 2021  − April 25, 2021],
T10	[460,491(31)]	[April 25, 2021  − May 26, 2021],
T11	[491,515(24)]	[May 26, 2021  − June 19, 2021],
T10∗	[460,515(55)]	[April 25, 2021  − June 19, 2021],
T11∗	[515,532(17)]	[June 19, 2021  − July 6, 2021],
T11∗∗	[515,627(112)]	[June 19, 2021  − October 9, 2021].

Fig. 3

(a) The rate v and its linear trend vs t in time subintervals from Table 1; (b) acceleration a and its linear trend X vs t, (c) phase-plane curve a vs v; (d) phase-plane rectangular curve (piecewise continuous acceleration) a vs linear trend of the rate from Eq. (1), together with curve a vs the average rate v over 7 days; calculated from the statistical data [1a].

Fig. 4

The same as in Fig. 3, but calculated from the statistical data [1b].

A powerful tool for studying dynamical systems and processes is the use of phase space (PS) [15]. In our case, it is formed by the function v(t) and its finite difference derivative a(t). The set of points M in R 2 with Cartesian coordinates v(t) = (N(t + 1) − N(t))/N , N = 106, and a(t) = (v(t + 1) − v(t)) will be referred to as the phase space of the COVID-19 pandemic. The points of this space adequately reflect all the characteristic features of the COVID-19 pandemic in the time period we are interested in.

In the classical mechanics, the theory of differential equations, and the theory of dynamical systems, phase trajectories are smooth manifolds. In the case of the COVID-19 pandemic, they are no longer smooth. Let us show that the phase trajectories of the COVID-19 pandemic resemble the trajectories of Brownian motion and, therefore, can be considered fractal curves. For this purpose, we present the PS of the COVID-19 pandemic according to data [1a] and [1b] in Figs.3c and 4c, respectively.

The fractal dimensions calculated by us for the COVID-19 pandemic from data [1a] and [1b] are D = 1.4918 and D = 1.6757, respectively. These estimates indicate a high degree of randomness of the spread of the COVID-19 pandemic over the time period under study. Figs. 3c and 4c clearly demonstrate the nondeterministic chaos of the COVID-19 pandemic behavior, since there are numerous self-intersections of the phase trajectory of this process. For comparison, in Figs. 3d and 4d we show the phase-plane rectangular curves (piecewise continuous acceleration) a vs linear rate trend from Eq. (1), together with curves of a vs average rate v over 7 days.

Multifractal dynamics of COVID-19 pandemic

It is assumed that and the curve v(t) is multifractal. The trend line has the unit fractal dimension and v(t) has the fractal dimension D . In the MFD mathematical model, the function X(D) is determined by the equation [13], [14].

The MFD parameters B 0 and η in the time t interval of interest are considered constant. The dynamics of in this model is determined by the analytical dependence of the coefficient A on the fractal dimension D. The values of D satisfy the condition 1 < D < 2. We divide this interval into two intervals 1 < D < D 0 and D 0 < D < 2. The parameter η in (4) describes the effective influence of basic factors on the spreading of COVID-19 (social restrictions, isolation, lockdown, vaccination, etc.). In Eq. (4) it is assumed that ∣X(D) ∣ < < 1 and ∣η ∣ < < 1, which corresponds to a sufficiently slow (quasi-equilibrium) character of the process dynamics.

The quantities B 0 and η will be considered constant during the time of the process observation, and the process dynamics is mainly determined by the analytical dependence of the coefficient A on the fractal dimension D. We use the following analytical representation for X(D):In Eqs. (4), (5), the MFD model parameters D 0 and D satisfy the condition 1 < D 0 < D < 2. At the point D 0, the functions in Eq. (5) are smoothly matched.

Remark. We emphasize, that the MFD parameters D 0, D , and η are found from the condition of minimum deviation of the calculated values X(D), e.g. X (D ), from Eqs. (4), (5), (6), (7) with respect to the statistical data [1a] (or [1b]) [16], as X determined by Eq. (1) at different intervals of the curve v(t).

We are interested only in real-valued solutions of Eq. (4) A cubic equation has either one or three solutions. We reduce Eq. (4) to the form more convenient for analysis. After a substitution X = (η/B 0)1/3 ξ to Eq. (4) it takes the form

The dependence of the solutions ξ of Eq. (6) on the parameter λ is shown in Fig. 5 . From Fig. 5 and Eq. (6), it follows that the point λ = λ(D ) = (27/3)1/3 is a bifurcation point. At λ < λ Eq. (6) determines one function ξ 1(λ) and at λ > λ - three functions ξ 1,2.3(λ). From Eq. (6) it follows that functions ξ 1,2.3(λ) have asymptotes:The functions ξ 1,2.3(λ) determine the stationary states of the dynamical system. If λ > λ , there are three possible states. Jumps between these states can occur under the action of nonstationary factors (social restrictions, isolation, lockdown, vaccination, etc.) on the dynamical system or under the loss of stability of the state described by the function ξ 2(λ) and realized before the jump. According to the statistical data [1a], the jumps of v(t) are observed in the dynamics of COVID-19 pandemic.

Fig. 5

Dependence of the solutions ξ of Eq. (6) on the parameter λ.

Let us divide the time interval 1 ≤ t ≤ 491 into ten intervals T 1 = [1, 50], T 2 = [50, 75], T 3 = [75,125], T 4 = [125,165], T 5 = [165,255], T 6 = [255,285], T 7 = [285,355], T 8 = [355,390], T 9 = [390,460], T 10 = [460,491] and construct ten piecewise linear trends , i = 1, 2, …, 10, on them. The division of the observation interval into parts T and the calculation of the parameters of the functions were carried out from the condition of the best approximation of the function v(t) by a piecewise linear function in the space L 2. For these intervals, the values of the parameter X of the linear trend rate , the fractal dimension D , and the parameter λ , calculated in the MFD model from Eqs. (4), (5), (6), (7) using the MFD algorithm presented below, are presented in Table 2 for COVID-19 data [1a] and [1b], respectively.

Table 2

The values of the parameter X - the linear trend of the rate , the fractal dimension D and the parameter λ, calculated in the MFD model from Eq. (4) for COVID-19 data [1a] and [1b]. Notation: (x) = 10.

[1a]
i	Ti (days)	Xi	Di	λi
1	1−50	0.3451(−4)	1.1460	32.9702
2	50−75	0.3418(−2)	1.0709	9.0274
3	75−125	0.3440(−3)	1.3359	3.3076
4	125−165	0.2369(−2)	1.2562	4.3344
5	165−255	0.9144(−3)	1.4119	116.4672
6	255−285	0.7704(−2)	1.1873	13.8228
7	285−355	0.1762(−2)	1.2135	60.4490
8	355−390	−0.9772(−2)	1.1659	7.4877
9	390−460	0.7643(−2)	1.3135	−9.5730
10	460−491	−0.1285(−1)	1.1236	7.7561

In the region of jumps, the piecewise linear trend coincides with function v(t). The plots of functions v(t) and are presented in Figs. 3a and 4a.

For the jumps of v(t), the value of X in the region of the jump X can be estimated by the formulawhere v is the value of v(t) in the maximum of the jump, v 0 is the value v(t) before the jump, τ is the jump duration. In the case of jumps of v(t) for the data [1a] shown in Fig. 3a we find X = 0.006533, X = 0.8281, X = − 0.2755. From (7) at λ 1,7,10 > > 1 we get the relations

Using our programs for the MFD model, implemented in Maple, we have shown that the curves , i = 1, 2, …, 10 are close to fractals with an accuracy not worse than 1 %. The values of fractal dimensions D for the ten considered intervals of time t are presented in Table 2 for the [1a] data.

For data [1a], let us consider the time intervals T = T 1 ∪ T 2 ∪ T 3, T = T 4, T = T 5 ∪ T 6 ∪ T 7, T = T 8 ∪ T 9 ∪ T 10. We begin the calculation of the MFD parameters D 0, D , η, B 0 from the time interval T . From Table 2a the inequalities followand from the formula (5) and the condition 1 < D 0 < D < 2 it follows that

Then the following system of equations can be written to determine the MFD parameters D 0, D , η, B 0:

The analytical solutions of the system of Eq. (12) have the form:

From Eq. (13) we get the values of the MFD parameters in the time interval T :

From Eq. (15), it follows that the condition for a big jump on December 10, 2020 (v = (70585238 − 69086805) = 1499433 = 0.15 ∗ 106) in the seventh time interval T 7 is fulfilled: λ 7/λ = 31, 9856 > > 1. In a similar way, the calculations were performed in the time intervals T , T , T . The resulting values of the parameter λ , i = 1, 2, …, 10, are presented in Table 2a, and the MFD parameters for periods T , T , T , T - in Table 3a for data [1a]. The calculations were carried out similarly in the time intervals T , T , T , T for data [1b]; the results are presented in Tables 2b and 3b.

Table 3

Parameters X, D0, D, η, B0 of MFD model for the COVID-19 data [1a] and [1b] for periods T, T, T, T: T = T1 ∪ T2 ∪ T3, T = T4, T = T5 ∪ T6 ∪ T7, T = T8 ∪ T9 ∪ T10. Notation: (x) = 10.

[1a]
	TI	TII	TIII	TIV
Xp	0.6533(−2)	0.8280	−0.2776
D0	1.1743	1.1743	1.2214	1.3002
Dk	1.2770	1.2336	1.2253	1.3154
η	0.1219(−2)	0.1219(−2)	0.2243	−0.7273(−1)
B0	−8.2750(5)	−8.2750(5)	−185.6646	−73.4681

The MFD mathematical model, as follows from Eq. (13), accurately determines the parameters of the trend and the big jump of the function v(t) for data [1a], describing the dynamics of the COVID-19 pandemic from day 1 to day 491 from its beginning (Fig. 2).

It follows from the data of Table 2 that the values of the parameter X - the rate of the linear trend v (t) for all time periods T , except the eighth and tenth, are positive. This means a monotonous growth of function v (t) until January 10, 2021, after which there was a qualitative change in the dynamics of the COVID-19 pandemic: periods of decline alternate with periods of growth. Thus, during the eighth time period of 35 days (from 10 January 2021 till 14 February 2021), the number of daily COVID-19 diseases in the world decreased from 0.7969 to 0.3006 million people per day. During the ninth period of 70 days(from 14 February 2021 till 25 April 2021), it increased by about the same value. For the tenth time period T 10 ∗ of 56 days (April 25, 2021 - June 19, 2021) a decrease in the number of daily infected people occurred again.

If the above pattern persists, then the expected forecast could be as follows: the eleventh time period T 11 (with both increasing and decreasing patterns) will last 112 days and during it, the function v 11(t) should also increase from about 0.3 to 0.8 million people daily. The subsequent periods are expected to demonstrate an oscillatory behavior. This forecast of the third wave based on the analysis of our calculations performed in July 2021 using the statistical data [1a] and [1b] is shown in Figs. 2 (left), Fig. 3, Fig. 4. It is confirmed by the current statistical data [1bc] [17] shown in Fig. 2 (right), which was added to the final version of the paper prepared in October 2021. Note, that based on the approximation of v(t) by the product of a power function and sint/T function (where T is the oscillation period), we get a forecast for October 10, 2021 v = 0.30 ± 0.11 million people, which is in good agreement with the data of daily statistics for this date, v = 0.30, and for October 11, 2021, v = 0.38 [17].

Let us supplement the range of our studies with a time interval from May 27, 2021 to December 26, 2021. We approximate the data of [16(a)] on this time interval with a piecewise linear trend consisting of three intervals: T 11=[512,567(55)] ([June 16, 2021 - August 10, 2021]); T 12=[568,629(61)] ([Aug 11, 2021 - Oct 11, 2021]); T 13=[630,705(75)] ([Oct 12, 2021 - Dec 26, 2021]). On these intervals, the parameters X , X and D turned out to be respectively equal: X 11 = 0.6164 ⋅ 10−2, X 12 = − 0.5086 ⋅ 10−2, X 13 = 0.4171 ⋅ 10−2, X = 0.31, D 11 = 1.1537, D 12 = 1.2316, D 13 = 1.3580. Using the relation D 11 < D 12 < D 13, solving a system of equations similar to systems (12) and (13) on the time interval T = T 11 ∪ T 12 ∪ T 13, we obtain the values of the MFD parameters on this interval: λ 11 = 11.9622, λ 12 = − 14.4961, λ 13 = 17.6773; D 0 = 1.2126, D = 1.2438, η = 0.1047. The inequality λ 13 > λ implies the possibility of large jumps in v(t) at t > T 13, which is confirmed by world statistics [16(a)].

Multifractal dynamics of basic reproduction number

The finite-difference equation for the number of infected people N(t) at time t in the reduced SIR model has the form [20], [21], [22].where N is the number of population, τ is the time of possible spreading the infection by a virus carrier, α(t) is the probability of infecting a healthy population member upon a contact with all population members per unit time (Δt = 1 day). This probability is defined as the probability of infection upon a single contact with a virus carrier multiplied by the number of his contacts with all members of the population per unit time. From Eq. (15), the dimensionless basic reproduction number α (t) = α(t)τ can be expressed as.and the rate of its daily change as Δα (t) = α (t + 1) − α (t):

Fig. 6 presents the basic reproduction number α (t) = α(t)τ for τ = 7 and the daily variation of the basic reproduction number Δα (t) = α (t + 1) − α (t) plotted from Eqs. (16), (17) using the statistical data [1a] and [1b] of the COVID-19 world dynamics.

Fig. 6

The basic reproduction number α(t) vs t from Eq. (16) (a, b) and its rate Δα(t) vs t from Eq. (17) (c, d) calculated using data [1a], [1b], respectively.

To estimate the parameters of multifractal dynamics and the fractal dimensions D in different time intervals T of variation of the basic reproduction number α (t), using the partitioning into a linear and rapidly oscillating trends we use the algorithm and programs from the previous Section 3. Table 4, Table 5 for the data [1a] and [1b] present the calculated parameters of α (t) multifractal dynamics shown in Fig. 6. The comparison of data in the tables shows that there is some difference in the structure of coefficients α (t) for the data [1a] and [1b]. Thus, according to Table 4, Table 5 [1a], there are three characteristic periods with different structure, while for the data [1b] there is only one such period. This is associated with the presence of three jumps of this coefficient for [1a] and only one jump for [1b]. In the time intervals of days 1 − 77, for the data of [1a] there are two jumps, and for the data of [1b] one jump. For the interval of days 77 − 491, there is one jump for [1a] and no jumps for [1b]. The structure of the parameter α (t) considerably differs for the cases [1a] and [1b]. In Table 5 [1a], the values of parameter η are zero in all time intervals, while in Table 5 [1b] they are nonzero in both time intervals. Worth noting is the substantial difference in the coefficients B 0 between two data sets by three orders of magnitude.

Table 4

Calculated parameters of the multifractal dynamics model for α(t) using data [1a] and [1b]. Notation: (x) = 10.

[1a]				[1b]
i	Ti (days)	Xi	Di	i	Ti (days)	Xi	Di
1	1–36	0.8626(−2)	1.1058	1	1–77	0.1100(−1)	1.1046
2	36–77	−0.9673(−2)	1.1242	2	77–491	0	1.6387
3	77–491	0	1.4252

Table 5

Calculated parameters of the multifractal dynamics model for α(t) using data [1a] and [1b].

[1a]				[1b]
i	1	2	3	i	1	2
Xp	5.0170	2.4519	1.3620	Xp	3.7584
D0	1.5600	1.5600	1.5600	D0	1.6387	1.6387
η	0	0	0	η	0.02059	0.02059
B0	−0.01805	−0.07248	−0.07248	B0	−12.0477	−12.0477

To estimate the daily increment of the basic reproduction number Δα (t) = α (t + 1) − α (t), we divide it into a linear trend and a fast-oscillating one, , and apply the algorithm [13] presented below. The sets of numerical values of the parameter Δα (t) are determined based on the statistical data [1a,1b] in the form of time series obtained from Eq. (17), which are approximated by piecewise linear curves. We cover these curves with cells of size h in R 2, use the program to find the number of cells M(h), and approximate M(h) with a power function of h.

The approximation parameters Γ and D are found from the condition of least deviation of functions M(h) and . Let us use the sequence h , k = 1, 2, …, K. Then the relative deviation δ of the set studied from a fractal is estimated using the formula

For particular calculations we took h = 1/(10k), K = 10. The results of our calculations of D and δ are presented in Table 6 for time periods T : T 1, T 2, T 3, T 4; t ∈ [t (,  t ] (days). From the data of Table 6 it follows that the daily increment of the basic reproduction number Δα (t) presented in Fig. 6 is a multifractal with the accuracy not worse than 0.95%, and the fractal dimensions D in different time intervals T vary within the range 1.0932 − 1.6033 for data [1a] and 1.1136 − 1.6183 for data [1b].

Table 6

MFD parameters D‘δ for Δα(t) in the intervals T1, T2, T3, T4 for data [1a] and [1b].

		[1a]		[1b]
i	Ti(days)	Di	δi	Di	δi
1	1–60	1.0932	0.003748	1.1136	0.003921
2	60–320	1.5837	0.005827	1.6183	0.005154
3	320–480	1.6033	0.004400	1.5601	0.009469
4	480–530	1.4245	0.008212	1.3357	0.005714

The dimensionless value of the basic reproduction number α (t) = α(t)τ at τ = 7, presented in Fig. 6, changes with time within the limits 0.7 − 2.5, which on average corresponds to the lower estimate 1.4 − 3.9 for COVID-19 and the corresponding lower estimate for the value of the probability rate α(t) is 0.1 − 0.35 (day−1).

In this case, there is an anticipatory change in its value α (t), presented in Fig. 6, in comparison with the change in v(t) = ΔN(t) = N(t + 1) − N(t) of the daily increase in the number of infected people at the moment of time t, shown in Fig. 2, Fig. 3, Fig. 4.

Conclusions

To describe the evolution of the function v(t), which represents the number of daily cases of COVID-19 infection in the world, we proposed using a mathematical model of multifractal dynamics. The model is free from hypotheses and assumptions about the structure of the process under study. This is an alternative to mathematical models based on differential and integro-differential equations. The multifractal character of the function v(t) is shown. The MFD parameters, were calculated and summarised in Table 2, Table 3, which precisely describe the piecewise linear trend, the great jump of function v(t) for data set [1a] and the absence of the jump for data set [1b]. Based on the analysis of our calculations performed in July 2021, a prognosis of the third wave of daily new cases of COVID-19 infection in the world. It is confirmed by the current statistical data shown in Fig. 2 (right) added to the final revised version of the paper, submitted for publication in October 2021.

Solving the finite-difference equation of the reduced SIR model with respect to the basic reproduction number ατ(t) and its increment Δατ(t) = ατ(t + 1) − ατ(t), we evaluated the fractal dimension parameters D for various intervals of their variation (these values of D lie within the interval from 1.1058 to 1.6387 for ατ(t) and within the interval from 1.0932 to 1.6183 for Δατ(t)) in the course of COVID-19 pandemic dynamics.

The proposed multifractal dynamics approach, resulting multifractal structures and interpretation provide additional clarity of vision and useful tools for short-term forecasting and identifying new trends in the dynamics of the COVID-19 pandemic, such as the new strain known as Omicron SARS–CoV-2. Note a starting number of infected Omicron SARS–CoV-2 doubled in three and a half days. What can one gives a forecast for the fourth pandemic wave? These data can be used as information about the beginning of the development of the fourth pandemic wave. Then Eq. (15) reduced to the following with a natural solution x(t) = x(t = 0) exp (αt). Comparing with x(t) = x(t = 0)2 one can obtain α = 0.198, i.e. 1000 infected people will infect 192 people in a day. The solution of (15) with α = 0.198 and τ = 14 gives a duration of the fourth pandemic wave equals 120–140 days, that is in agreement with COVID-19 World Statistics [16] (b).

Future applications of multifractal dynamics mathematical model to describing dynamic systems and processes, as well as the comparison of its forecasts with experiment, are expected to reveal the model viability.

CRediT authorship contribution statement

The authors declare that the study was realized in collaboration with the same responsibility. All authors read and approved the final manuscript.

Declaration of competing interest

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