Simulation of Spatial Spread of the COVID-19 Pandemic on the Basis of the Kinetic-Advection Model

Abstract

A new two-parameter kinetic equation model is proposed to describe the spatial spread of the virus in the current pandemic COVID-19. The migration of infection carriers from certain foci inherent in some countries is considered. The one-dimensional model is applied to three countries: Russia, Italy, and Chile. Both their geographical location and their particular shape stretching in the direction from the centers of infection (Moscow, Lombardy, and Santiago, respectively) make it possible to use such an approximation. The dynamic density of the infected is studied. Two parameters of the model are derived from known data. The first is the value of the average spreading rate associated with the transfer of infected persons in transport vehicles. The second is the frequency of the decrease in numbers of the infected as they move around the country, associated with the arrival of passengers at their destination. An analytical solution is obtained. Simple numerical methods are also used to perform a series of calculations. Calculations us to make some predictions, for example, about the time of recovery in Russia, if the beginning of recovery in Moscow is known.

1. Introduction

The COVID-19 pandemic is a serious challenge not only for medicine, but also for other scientific disciplines including physics and mathematics. Numerous scientific groups worldwide are working on this problem based on different and complex models using high performance computing.

There are many models that aim to describe epidemics. Let us note the well-known susceptible, infectious, and removed individuals (SIR) model and similar ones. In this sense, first, one can mention the original paper [1]. The vast review of a family of the methods is presented in the monograph [2]. Although the basic mechanisms of each phenomenon are different, their effective mathematical description often defines similar constitutive equations and dynamical behaviors framed in the general theory of reaction-diffusion processes [3], see also [4].

The family of SIR-based models describe the changes in the number of susceptible, infectious, and removed individuals over a considerably large period of time. Some more sophisticated models consider the specifics of COVID-19, although they become complicated due to the number of parameters that need calibrating [5]. For an overview see also [6]. Note that, in general, SIR models do not take into account the population density at different points in space. The main difference of our model from SIR-like models is that we intend to describe the spatial spread of the epidemic in time (which SIR models do not) when local infections can be neglected. We study the movement of the wave of infections along the selected direction, but we do not consider the entire time evolution of the epidemic in a given area, since a few weeks after the first infections, it is already determined by local infections.

Our model takes into account population density, which strongly determines the profile of the spatial infection wave. We should mention the model of the spatial propagation of the first wave of the virus in space and time [7]. As proposed in this model, the number of infected individuals is supposed to be proportional to the power of distance of the source of epidemics and the power of the population in a given city. The powers of the values are obtained using freely available data about cases of COVID-19. We also propose a two-parameter model, but these parameters have strict physical meaning: the speed of migration of infected individuals and resistance value that depends on various factors, for instance, the quarantine measures. We show that the exponential decrease of the number of infected individuals is a result of the transport equation, and it depends not only on distance but on the population density.

In [8], the investigation is presented how the number of cases by COVID-19 scale with the population of Brazilian cities. Here, e.g., urban scaling relations of COVID-19 cases, the time dependence of the scaling exponents for COVID-19 cases, and the association between growth rates and city size are shown. However, in that paper, the spatial spread of the epidemy over long distances is not studied. There is some statistical research on how different factors affect the COVID-19 spread (see [9,10,11]). The most relevant to this work is the analysis of spatial factors. However, in those papers, there is no proposition of a deterministic model that describes the spatial propagation of a virus around countries.

Various spatial models of the spread of viruses are known, but they often contain many parameters that require determining their exact mutual correspondence. Simpler models that convey qualitative and quantitative patterns of the distribution would certainly be of interest, both scientific and practical.

Among the factors contributing to the spread of the COVID-19 infection, one can note the current significant speed of vehicles and a large number of social contacts in cities. At the same time, modern information technologies make it possible to quickly track down the spread of the pandemic, as well as collect sufficiently accurate data to construct adequate mathematical models (the resources [12,13,14,15] are used). There are studies devoted to the problems of modeling biological and social processes. Note that in [16,17] a diffusion-type equation is used. On the basis of such methods, various autowave processes were studied, including the spread of infections, see, e.g., [18]. In [19], kinetic transport equations were considered to study the nature of traffic flows. In [20], kinetic equations were used to model some sociohistorical processes including a model of aggression.

It is important that in the dynamical changes of the spatial distribution of COVID-19, two independent factors coincide. One is the advection related to the transport vehicles, the other the counteraction with the outer media. In fact, the latter denotes decreasing the number of passengers coming back to their home. The adequate apparatus is the kinetic equation; the direct methods for this is described in [21]. However, in the present case the equation is quite simple; we suggest a one-dimensional advection-kinetic type model of the spatial movement of infected carriers in Russia, Italy, and Chile. For these countries, the migration of carriers of the virus occurred from the designated centers: Moscow (the spread of the virus eastward from the capital is confirmed), the province of Lombardy, in Italy, and Santiago, Chile, respectively. The use of a one-dimensional model is justified by socio-geographical conditions: the proportion of “length” to “width” is quite large. Lateral invasions of the virus can be ignored. In Russia, the northern areas are sparsely populated, in the south there are deserts and mountains on the borders with other countries. In Italy, from the west and east, the Apennine peninsula is outlined by the Tyrrhenian and Adriatic seas, respectively. In Chile, the Pacific Ocean is to the west, and the high Andes to the east.

The one-dimensional model uses averaging over rather large territories of countries. The map of each country is divided into large bands. For Russia, 7 bands are allocated, for Italy—9, for Chile—4 and 3 (to the north and south of Santiago, respectively). The proposed model is determined by the following parameters: $U$ is the average speed of transport coming from the center of infection and σ is the “resistance” to the movement of carriers of the virus by vehicles. In this kinetic-figurative approach we consider the early days of the COVID-19 pandemic in the spring of 2020. This model is close to the models used in [20]. But in the present work it is simpler since it is linear. In contrast to the mentioned work, where background elements that participated in the interaction were eliminated, in our model the background parameters do not change with the advancement of the infected elements in space.

Our model, in general, is a “ballistic-kinetic” rather than a “reaction-diffusion” model. The fundamental difference is that in our approach, the spread of infection is not associated with a gradient in the number of infected. In the known model, the infection moves from an area with a large number of cases to an area with a smaller number, in other words, according to the corresponding diffusion equation. In our method, we consider the spread, which is not directly related to the number of infections, since the spread rate is due to another factor, namely, the average speed of transport (including cars, buses, trains, and airplanes). Thus, the spread of infection can occur from a region with fewer infected people to a region with a large number of infected people, that is, in the direction of the infection gradient. In the situations under consideration, indeed, the initial concentration of diseases, for example, in Moscow, is much higher than in the periphery, but the spreading process cannot be described by diffusion.

Two mechanisms of the infection spread can be defined: a transferential mechanism and a contact one. In our work, the transferential mechanism is studied, that is, infections caused by the migration of infected carriers from the center of infection. The contact mechanism comes in after the relocation of infected carriers from the center, and it is due to the contagion within the regions. The superposition of these two factors gives the sum of all infected. The first phase of the process involves transferential contamination. This phase lays the foundations for the next phase of the infection development associated with local contacts. Contact infections are more pronounced with the introduction of quarantine isolation measures in individual regions. The third phase is associated with the spread of “recovery”, which is determined, in particular, by reaching a maximum of new infections per day for each region. Infected individuals travel by various vehicles. We consider the average speed, which takes into account the speeds of different air-, land-, and water-based vehicles. Another parameter of the suggested model is the value of “resistance” to the movement of virus carriers. The main equation of the model is the inhomogeneous transport equation with a kinetic term on the right-hand side: the decrease in number of carriers is taken into account which is associated with the carriers returning to their places of residence from the center of infection. Quarantine barriers between regions also play a role. Some predictions according to the present model are discussed, namely, the data of the recovery for the whole country, as we know that the velocity of recovery is approximately equal to the velocity of infection.

2. Mathematical Model

Analogously to [20], in order to determine the density of infected passengers disembarking at a certain place at a certain point in time, the following non-uniform transfer equation is applied:

$$\frac{\partial n}{\partial t}+U\frac{\partial n}{\partial x}=-\sigma n\left(t,x\right),$$

(1)

where t is the time, x is the distance, n(t,x) is the density of moving virus carriers in a vehicle, $U$ is the speed of a vehicle, and σ—the coefficient of “resistance” to the movement of infected elements mainly due to the landing of passengers in places of residence, it is measured in units of temporal frequency. This is a simple form of the well-known kinetic equation, in which the left-hand side describes the advection of particles, and on the right-hand side the collision process is presented, here we consider the linear case without the gain term. We emphasize that this is not an ordinary diffusion-reaction equation. Here, advection is not the result of a corresponding mass gradient as in the diffusion equation. The speed of movement of passengers (or elements or particles in physical terms) is related to the speed of the vehicle.

The initial condition for the Cauchy problem is:

$$n_0\left(x\right)=H\left(-x\right),$$

where H(x) is the Heaviside function. In other words, $n_0$ = 0 for t = 0 and for x > 0. This formulation of the problem means that in the studied area where x > 0 carriers of the virus enter across the border x = 0, which corresponds to the idea of the main permanent source of infection.

Let us use the method of characteristics to solve Equation (1). Let $x=x\left(s\right), t=t\left(s\right),$ thus:

$$\frac{d}{ds}n\left(t\left(s\right), x\left(s\right)\right)=\frac{\partial n}{\partial t}·\frac{dt}{ds}+\frac{\partial n}{\partial x}·\frac{dx}{ds}.$$

Comparing this expression with the left-hand part of Equation (1), we see that if $dt/ds=1$, $dx/ds=U$, we get the following ordinary differential equation:

$$\frac{dn}{ds}=-\sigma n\left(s\right),$$

which has an exact solution: $n=A{e}^{-\sigma s}$. From $dt/ds =1$ and $t\left(0\right)=0$, we get that $s=t$, thus from $dx/ds=U$, we obtain that $x=Ut+x_0$, where $x_0=\mathrm{const}$. Thus we get:

$$n\left(t,x\right)=A_1{e}^{-\sigma \frac{x}{U}}.$$

In fact, we use standard methods for solving linear partial differential equations, see, e.g., [22]. Assuming n(t,0) = 1 for t > 0, the solution is as follows:

$$n\left(t,x\right)=n_0\left(x-Ut\right)\cdot e^{-\sigma \frac{x}{U}}.$$

(2)

Here, a day is taken as a unit of time, a km is taken as a unit of distance, respectively, the speed is expressed as km/day, “resistance coefficient” (frequency) has the dimension ${\left(\mathrm{day}\right)}^{-1}$.

We will denote the density of virus carriers that landed at a given point by $n_M$(t,x). It is clear that this density increases with the decrease of n(t,x), thus we can write:

$$\frac{dn}{dt}+\frac{d{n}_M}{dt}=0,$$

and, since

$$\frac{dn}{dt}=\frac{\partial n}{\partial t}+U\frac{\partial n}{\partial x},$$

taking into account Equation (1), the equation for $n_M$(t,x) is:

$$\frac{\partial n_M}{\partial t}=\sigma n.$$

(3)

The solution of this equation is given by:

$$n_M\left(t,x\right)=n_{M_0}\left(x\right)+\underset{0}{\overset{t}{{\displaystyle \int }}}\sigma ·n\left(\tau ,x\right)d\tau .$$

The substitution of Expression (2) and the zero initial condition allows us to an expression for $n_M$(t,x):

$$n_M\left(t,x\right)=\frac{\sigma }{U}·n_0\left(Ut-x\right)·e^{-\sigma \frac{x}{U}}.$$

(4)

The density of disembarked passengers $n_M \left(t,x\right)$ grows with time as infected new people arrive. It is important to note that the profile of infected passengers in the vehicle n(t, x) with time goes to some definite curve, exponentially decreasing in space. However, the density $n_M\left(t,x\right)$ will constantly increase in time at each spatial point if conditions are not changed. The introduction of an area isolation regime can be an example of such change. In this case, it is necessary to change the boundary conditions.

We also used a simple numerical method to obtain a large number of variants of a solution. We applied the finite difference scheme [23].

In general, σ can be a function of t and x. $U$ can also be a dependent quantity, but we assumed it to be constant. One of the options for the dependence of σ on the coordinate can be obtained by taking into account the population density:

$$\sigma \left(x\right)={\sigma }_1·\rho \left(x\right),$$

(5)

where ρ(x) is the linear population density, having dimension 1/km, and ${\sigma }_1$ is an arbitrary constant of dimension km/day.

Substituting Equation (5) into Equation (2), we obtain the following expression:

$$n\left(t,x\right)=n_0\left(x-Ut\right)·e^{-\frac{{\sigma }_1}{U}{{\displaystyle \int }}_0^x\rho \left(z\right)dz}.$$

A similar substitution into Equation (4) yields

$$n_M\left(t,x\right) = \frac{{\sigma }_1·\rho \left(x\right)}{U}·n_0\left(Ut-x\right)·e^{-\frac{{\sigma }_1}{U}{{\displaystyle \int }}_0^x\rho \left(z\right)dz}.$$

(6)

In order to determine the parameter ${\sigma }_1$, we fix t in Equation (6) and write it down for two values of x: x = 0, the center of infection spread and x = Δ, the distance to the second band, and approximate the integral by the trapezoid formula:

$${\sigma }_1=\frac{2U}{\left({\rho }_1+{\rho }_2\right)\mathsf{\Delta }}\ln \frac{{\rho }_2{n}_M\left(t,0\right)}{{\rho }_1{n}_M\left(t, \mathsf{\Delta }\right)}.$$

(7)

As noted above, along the selected direction of the spread of infection, the territory of the country is divided into stripes. In open sources there is information on the number of new cases of infection and recovery for each day by region [12,13,14,15]. To approximate these data in the selected bands, we summarized the values in the regions included in the band with the corresponding weights. The weight was calculated as the ratio of the population density in the region to the sum of the densities for all the regions included in the strip. Let $a$ be the number of new infections per day for a strip. For the i-th region, the number of new infections $a_i$ is known. Suppose that the strip includes N regions. We have a set of values $a_1$,…,$a_n$. Len ${\rho }_1$,…,${\rho }_n$ be the population density in the regions included in the strip. Then the parameter a can be represented as a weighted sum:

$$a=\frac{{{\displaystyle \sum }}_{i=1}^N{\rho }_i{a}_i}{{{\displaystyle \sum }}_{i=1}^N{\rho }_i}.$$

Russia, Italy, and Chile were selected to apply the model. In these countries, the spread of infection occurred mainly in one direction from a dedicated center. This was facilitated by the uniqueness of the situation and the geographical characteristics of the countries. The most important sources of the epidemic were associated with the metropolitan centers of the countries (Moscow—the capital of Russia, Milan—the capital of Lombardy in Italy, Santiago—the capital of Chile), which were visited by tourists or citizens of this country, returning from other countries in which the epidemic had already begun. In Russia, citizens returning from abroad go through Moscow to the east to their places of residence. The geographical extent of the country allows us to consider a one-dimensional equation.

3. Simulation of the Spread of the Epidemic in Italy

In Italy, there was a spread from the Lombardy region to the southeast along the narrow Apennine peninsula. The country’s territory was divided into 9 bands (

Table 1. Regions to bands correspondence in Italy.

№	Distance, km	Regions
1	0	Lombardy
2	166	Piedmont, Veneto, Emilia-Romagna, Liguria, Aosta Valley, Trentino-South Tyrol
3	272	Friuli Venezia Giulia, Tuscany
4	376	Marche, Umbria
5	482	Abruzzo, Lazio
6	592	Molise
7	690	Campania
8	781	Apulia, Basilicata
9	933	Calabria

), in which the values for a number of infected people were found as weighted sums for the regions included in the band. The population density is shown in Figure 1 where the distribution in the bands is taken into account.

The $U$ parameter, which determines the rate of spread of infected elements, can be estimated from the data on the first infections in each band. We used the first day, on which there were at least 5 cases. Having constructed the points corresponding to these days, using the least squares method, we will construct a straight line passing through the point corresponding to the center of infection (Lombardy). $U$ is equal to the tangent of the slope of this straight line in Figure 2a.

The $U$ parameter, which determines the rate of spread of infected elements, can be estimated from the data on the first infections in each band. We used the first day, on which there were at least 5 cases. 5 is an arbitrary number, but the higher the number the more it is affected by the contact mechanism of the virus spread, and the lower the number, the more it is affected by the statistical error. Having constructed the points corresponding to these days, using the least squares method, we will construct a straight line passing through the point corresponding to the center of infection (Lombardy). The additional condition is that the line passes through the origin point which corresponds to the day when there were at least 5 cases in the first band. $U$ is equal to the tangent of the slope of this straight line in Figure 2a.

Although the average speed of vehicles can be obtained through other methods or by using a higher order polynomial for the least squares method, the method described here, despite its simplicity, can yield some prediction power which is shown in the results section of this paper.

The values from Figure 2b were used to determine the speed of recovery. The way it was done is explained later in the text.

For Italy, we get the speed: $U$ = 53 km/d. Using Equation (7), we obtain ${\sigma }_1$ = 0.0012 km/day.

A comparison of the numerical solution for $n_M$ with real data for Italy at the initial stage of infection is shown in Figure 3. The graphs show that the nature of the real and theoretical distributions is similar to each other, but the values of the real data are higher than the theoretical ones. This phenomenon is observed due to the larger population in the second band, which is due to the large number of Italian regions included in the second band.

Here is a consideration when comparing population density and $n_M$ distribution in Italy, see Figure 1 and Figure 3. There is a certain correlation between extrema on the two graphs. Moreover, the second maximum in Figure 3 is small due to a decrease in the number of moving infected. It is important that these maximums on the real distribution graph and the calculation graph correspond approximately to each other. There is a characteristic detail in both graphs, namely, the data maxima are shifted to the left compared to the maxima on the population density graph, which can be explained by the peculiarity of the decrease in the number of infected passengers. It can be assumed that the maximum of infected is associated with the maximum of the derivative of the density. The increase in population density is associated with the appearance of cities in these places on the map and, consequently, with the train stations, bus stations, and airports on the border of these cities. On the other hand, quarantine barriers can also play a role.

We define the onset of recovery for the selected band as reaching the maximum infections per day. The least squares method was used with an approximation using a second order polynomial. For each parabola, a maximum was determined, which correlated with the date of the beginning of recovery for this band. The data on the maximum infections per day for each band is shown in Figure 4 (only some bands are shown for brevity). In these plots, 0 on the x-axis is 1 March.

The rate of recovery was determined in the same way as the rate of infection (using the method of least squares with a linear approximation), see Figure 2b, the $U$ amounted to 50 km/d. According to the schedule of the spread of recovery, one can predict the approximate time of the beginning of recovery for Italy as a whole, namely the date when the line of distribution of recovery passes the point corresponding to about half of the population of the country. In Figure 2b, the origin on the abscissa corresponds to the beginning of recovery (reaching the maximum infections) in Lombardy. Due to the heterogeneity of the density distribution (see Figure 1), approximately half of the population of Italy corresponds to a territory up to a distance of 600 km from Lombardy. This point in Figure 2b corresponds to 2 April. As predicted from Figure 2b, the time of the beginning of recovery for the whole of Italy is corresponding to the beginning of recovery according to Figure 4, which shows that the maximums of infection are reached in the first week of April.

The rates of infection and recovery were approximately the same. This conclusion is consistent with the assumption that the conditions for the course of the disease were akin to each other in different areas, which is ensured by strict quarantine isolation of each area.

4. Simulation of the Spread of the Epidemic in Chile

In Chile, the main source of the spread of the virus was the capital of the country, Santiago, located approximately in the middle of the extended country. Four bands were allocated north of Santiago and 3 bands south of the capital. The data for Chile were aver-aged for each band by the regions included in it the same way as for Italy.

In the simulation, two different directions are considered independent of each other, that is, it is assumed that the number of infected passengers who left, for example, in the direction of the north, and then turned south, will be negligible compared to the total number of passengers. The bands in the north direction are shown in

Table 2. Regions to bands correspondence to the north of Santiago in Chile.

№	Distance, km	Regions
0	0	Santiago
1	50	Coquimbo, Valparaíso
2	470	Atacama
3	850	Antofagasta
4	1340	Tarapacá, Arica and Parinacota

. The bands in the south direction are shown in

Table 3. Regions to bands correspondence to the south of Santiago in Chile.

№	Distance, km	Regions
0	0	Santiago
1	50	O’Higgins, Maule, Ñuble
2	470	Biobío, Araucanía, Los Ríos
3	850	Los Lagos

.

Further to the south of Los Lagos, there are two more regions. These regions are not taken into account in this work, because of low population density in these regions so the statistics on coronavirus are comparable to the error in other regions.

The population density is shown in Figure 5. In the graphs hereinafter, Santiago is not considered (reference point x = 0), since the population density in this city is an order of magnitude higher than in other bands. Therefore, the density is taken into account from the first band in each direction—they are located at a distance of 50 km to the north and south of Santiago. From Figure 5 it can be seen that to the north, the population density from the first strip begins to decrease and then rises, and in the south, the density initially increases slightly, and then decreases.

The parameter values obtained for Chile by the same methods as for Italy have the following values: $U$ = 44 km/d, ${\sigma }_1$ = 0.0026 km/day for the north and $U$ = 97 km/d, ${\sigma }_1$ = 0.0065 km/day for the south. The appropriate graphs for southern Chile are shown in Figure 6 (for the north they are similar).

A comparison of the numerical solution for $n_M$ with real data at the initial stage of infection for southern Chile is shown in Figure 7. Figure 7a,b have some similarities, despite the rather uneven growth of real data.

The rates of recovery were found the same way as for Italy (Figure 8). The rates of recovery: $U$ = 23 km/d to the north of Santiago and $U$ = 96 km/d to the south. A mismatch in the rate of recovery of the rate of infection in the north may indicate some inaccuracy in the data or data insufficiency. Again, in Figure 8 only data for the south is shown and only plots for one band and all bands to the south of Santiago are shown for brevity.

In Figure 6b the abscissa origin corresponds to the beginning of recovery (reaching the maximum infection) in Santiago. Despite different distributions of population density to the north and to the south of Santiago, about half of the population in each direction corresponds to a distance of 400 km from the capital (see Figure 5). This point in Figure 6b corresponds to July 7 in the south. Predicted from Figure 6b, the prediction for the south turned out to be less accurate than it was for Italy.

Plotting for the north a graph similar to Figure 6b and applying the same method, we would get 12 July. Plotting for north a graph similar to Figure 8b would show that the maximum of infections per day is reached in mid-July which means a good correspondence between two different methods.

5. Simulation of the Spread of the Epidemic in Russia

In the direction to the east, the territory of Russia on the map was divided into 7 bands, as given in

Table 4. Regions to bands correspondence in Russia.

№	Distance, km	Regions
	0	Moscow, Moscow region
2	500	Nizhny Novgorod Region, Ivanovo Region, Republic of Mordovia, Penza Region, Tambov Region, Volgograd Region, Ulyanovsk Region, Kostroma Region, Republic of Chuvashia, Republic of Mari El, Saratov Region
3	1160	Perm Territory, Republic of Komi, Republic of Bashkortostan, Orenburg Region, Samara Region, Republic of Tatarstan, Kirov Region, Nenets Autonomous District
4	1850	Tyumen region, Khanty-Mansi Autonomous Okrug, Kurgan region, Sverdlovsk region, Chelyabinsk region
5	2150	Omsk Region, Tomsk Region, Yamalo-Nenets Autonomous District, Novosibirsk Region
6	3600	Krasnoyarsk Territory, Irkutsk Region
7	6000	Khabarovsk Territory, Jewish Autonomous Region, Amur Re-gion, Primorsky Territory

. The population density by bands is shown in Figure 9a. The source code for computing the weighted values of infected individuals and for computing $n_M$ is available in Appendix A.

The parameters obtained for Russia by the same methods as for Italy have the following values: $U$ = 87 km/d, ${\sigma }_1$ = $0.0031$ km/day. The graph according to which the speed was found is shown in Figure 10a. The speed turned out to be higher than the same parameter in Italy due to a much greater length of the country, hence the average speed of vehicles is higher. For example, it can be assumed that airplanes (the fastest mode of transport) are used for domestic traffic in Russia more often than in Italy, which as a result gives a higher travel speed in Russia.

Numerical calculations for the time after two weeks of infection confirm qualitatively the development of the infection, and, naturally, lower values compared to the real data. Here the second—the contact factor—plays the main role. It is assumed that the quarantine conditions are strictly fulfilled, so the transfer of infection does not matter anymore.

In Figure 9b, it can be seen that the real values for two consecutive dates begin to grow faster than the calculated ones due to infections in the field.

The rate of recovery was found the same way as for Italy and amounted to 52 km/day.

In Figure 10b, the x origin corresponds to the beginning of recovery (reaching the maximum infection rate) in Moscow. About half of the Russian population corresponds to a distance of 1000 km from Moscow (see Figure 9a). This point in Figure 10b corresponds to May 30. As predicted from Figure 10b, the time of the beginning of recovery for the whole of Russia is corresponding to the beginning of recovery according to Figure 11b, which shows that the maximum of daily infected reached by the beginning of June.

Currently, a new wave of the pandemic has emerged, and such clear conditions for the development and spread of COVID-19 are no longer observed. The autumn wave of the pandemic is more complex, which was facilitated by many factors, in particular, there was no such strict isolation as in the spring. Nevertheless, it can be assumed that Moscow was the main center for the spread of the virus to the east. Therefore, the maximum infection in Moscow can be reached about three weeks earlier than in the whole of Russia (see Figure 12). Our prediction concerning the difference between the mentioned maximum in Moscow and in Russia as whole: in Moscow on around 10 December and we expect that the maximum in Russia will be around 31 December.

In Figure 13 some additional data known in January 2021 was added to check the prediction. The maximums of parabolas moved by about 2 weeks to earlier dates (this is a consequence of the method of the least square root) but the distance between maximum in Moscow and maximum in Russia stayed approximately the same.

6. Conclusions

This paper presents a model capable of replicating the spread of coronavirus (COVID-19) in some countries with a distinct focus on the disease. A two-parameter kinetic one-dimensional model of the spatial spread of coronavirus infection was proposed using the example of three countries: Russia, Italy, and Chile. An essential assumption is related to introducing strict quarantine starting from a certain time and fulfilling it for a sufficiently long period. It was assumed that, before the start of quarantine, the transmission of infection took place practically only due to returning passengers from the mentioned centers, the way it happened in three countries under consideration. A comparison of the calculations with real data showed good agreement for the first weeks of the transmission of infection.

The process of recovery is in a sense “symmetrical” to the process of infection: recovery begins earlier in those places where the infection occurred earlier. Thus, the results obtained also have predictive capabilities. By contrast, the second wave of the pandemic does not fit into such a clear picture. The autumn wave is more complex, which was facilitated by many factors (in particular, there was no such strict isolation as in spring 2020). Nevertheless, it can be assumed that, during the second wave, Moscow was again the main center for the spread of the virus to the east of Russia. Therefore, the maximum infection in Moscow can be reached about two or three weeks earlier than in the whole of Russia. Our important forecast points to a difference between the indicated maximum in Moscow and Russia as a whole: in Moscow on around December 10, and we expect the maximum in Russia to be closer to the end of December. This prediction was confirmed taking into account new data from January.

Further work can be done with this model or its modifications. For example, it is possible to consider the migration of virus carriers returning through Khabarovsk, which is located, on the Far East of Russia, in addition to Moscow. This would put additional boundary conditions and thus alter the resulting equation. A similar modification of the problem for Italy would make sense not on the right (southeastern) border, but in the middle of the country for Rome and Naples.

The study of the more detail spread of the epidemy needs using “2D maps”, i.e., using two-dimensional form of the equation. From Moscow the waves of the epidemic spread along radial rays. In the simplest case of spreading the wave from Moscow to the east one can consider a central symmetry, the radius of the circle will be the only independent variable.

Strictly speaking, there was the other important center of infection in Russia, namely, in St. Petersburg, to the west of Moscow, but the spread of epidemics westward from the capital was not examined. The area between Saint Petersburg and Moscow can be considered as a two-dimensional region with an advection-kinetic equation without any symmetry. The model in a two-dimensional form with many centers of infection requires more complex numerical calculations (it is a significant fact that the primary source of infection appeared in St. Petersburg earlier than in Moscow).