Particle modeling of the spreading of coronavirus disease (COVID-19)

By the end of July 2020, the COVID-19 pandemic had infected more than 17 × 106 people and had spread to almost all countries worldwide. In response, many countries all over the world have used different methods to reduce the infection rate, such as case isolation, closure of schools and universities, banning public events, and forcing social distancing, including local and national lockdowns. In our work, we use a Monte Carlo based algorithm to predict the virus infection rate for different population densities using the most recent epidemic data. We test the spread of the coronavirus using three different lockdown models and eight various combinations of constraints, which allow us to examine the efficiency of each model and constraint. In this paper, we have tested three different time-cyclic patterns of no-restriction/lockdown patterns. This model’s main prediction is that a cyclic schedule of no-restrictions/lockdowns that contains at least ten days of lockdown for each time cycle can help control the virus infection. In particular, this model reduces the infection rate when accompanied by social distancing and complete isolation of symptomatic patients.


I. INTRODUCTION
Statistical mechanics provides a set of very powerful tools to model various biological and medical problems (see, for example, Refs. 1-3 and many more). One of the most studied these days is diffusion of pandemics, prompted by the current COVID-19 emergency (see, for example, Ref. 4). In addition, much of current research studies the physical aspects of the spreading of the virus (see, for example, Refs. [5][6][7][8]. Many techniques currently employed are based on the solution of differential equations (see, for example, Refs. 9-11 and many more) or fitting formulas (see, for example, Refs. 12 and 13). Both techniques are based on varied parameters to obtain several scenarios that are then treated as parts of a statistical ensemble for analysis. For example, in many countries (e.g., Germany and Italy), there is ample discussion about the role of the so-called R 0 parameter, i.e., the average number of individuals that a single actively infectious person can pass the virus to. The procedures to estimate R 0 are all based on a posteriori analyses but are usually part of the parameters that governments use to decide on measures to be taken.
In this paper, we propose a method that is essentially based on modeling a population as a set of interacting classical particles, each one with three states relative to the health status (susceptible to infection, infected and contagious, and recovered/died), in which standard thermodynamical parameters (temperature and density) are used to describe the characteristics of the population. This method allows us to apply the model to very different situations, ranging from a city suburb population to a single university classroom. The algorithm is based on standard Monte Carlo (MC) procedures of sampling the transition among subsequent states, which are essentially sampled from a statistical distribution, in the spirit of transport MC algorithms.

ARTICLE scitation.org/journal/phf
In our model, a healthy person (i) can become sick with a daily probability, Pi = ∑jPij, where Pij is a function of the distance between each infected person (j) in the area and the healthy person (i). We are treating the coronavirus spread as a "one-way" Monte Carlo Ising model as follows: A healthy person becomes sick as a result of an interaction with a sick person (or people), but a sick person stops being sick (i.e., recovers or dies) within an average time of ∼14 days although up to 40 days for severe cases (see Ref. 14 for the epidemiology data). After that time, the recovered person can no longer infect another person and cannot be sick again.
In contrast to other infection models, such as SIR [9][10][11] in this approach, the parameter R 0 is a direct outcome of the simulation and not pre-assumed. This is achieved by the relation between R 0 and the doubling time T d , which is a direct result of the infection probability chosen, which is, in turn, a function of observable epidemiological data and features of the studied population (e.g., average density on a given area or mobility). In the following, we will present the results of several simulations meant to reproduce the spread of the coronavirus in the presence of different lockdown constraints. The model's high flexibility enables us to control many parameters such as social distancing (SD), infection from an unknown source, etc. In Sec. II, we will describe some details of the model. In Sec. III, the different models of lockdown considered will be discussed. Section IV is devoted to the presentation and discussion of the results, and Sec. V, to the conclusions.

II. THE PARAMETERS AND PRELIMINARY ASSUMPTIONS
Given the scarce information available, we had to make some assumptions based on current data, which may be more stringent than it might be required by the real nature of the virus. In particular, we relied on Ref. 14 for the coronavirus epidemiology data and on Refs. 5 and 6 for the physical properties of the virus. We are aware that this model cannot take into account every single spreading event, but since such events affect the initial T d (which is a function of the population density) and since the infection process is random, we expect that the existence of such events will be reflected in the numerical results. In addition, in contrast to real life, here, there is no time gap between getting infected and being tested positive for the coronavirus. Therefore, an immediate decrease in the rate of infection resulting from lockdown is expected in the model, in contrast to the real data (see Ref. 15 for up-to-date data).
All simulations are performed assuming a surface area unit of 1 km 2 . Periodic boundary conditions are used, allowing us to get rid of broad confinement effects (e.g., the lockdown of an entire province or city) and look at the local dynamics of the infections within that area. In this work, we modeled the spread of the COVID-19 virus as a function of the population density in a specific surface. The application of periodic boundary conditions means that we have an infinite number of identical systems; each system is a replica of the others, i.e., if a person leaves the simulation surface on one side, an identical person will enter the surface from the other side. Population density is a function of the number of households in a certain area since it is crucial to distinguish between the infection among household and non-household contacts. 16

A. Parameterization of the model
We have identified a list of parameters and corresponding values that describe the population and the infection's kinetics. This list is obviously partial, but it could be quite easily extended.
1. The probability of developing symptoms over time t. This is described by a Gaussian peaked att = 5 days and with a standard deviation σt = 1 day. 2. The number of effective households, denoted by N. 3. The fraction of "silent carriers," which have no symptoms (aka asymptomatic) but can infect other people. Their fraction in the population is denoted by a silent , and the probability of transmitting the infection has been set to 0.5. 4. Each sick person is considered contagious between the third and the seventh day.
Simulations are started with a single infected person (the zero patient). In some runs, infections from an unknown source (a healthy person becomes sick without interaction with a known sick person) are allowed. For the first 14 days, the infection has not been detected yet, and the population walks freely without any restrictions. In some of the simulations, we "force" sick people with symptoms to maintain a distance of 8 m (i.e., stay at home) after day 14. This restriction reduces the probability of non-household infection.

B. Population dynamics
Our model is based on the principles of Brownian motion such that for each day, the population position (R) and displacement (ΔR) are given by where ΔR = √ Δx 2 + Δy 2 is distributed normally, where Δx (Δy) is the displacement in the x-(y-) direction and σ 2 R , the variance, is a function of the diffusion constant, D, where t = 1 day. For a Brownian motion, the diffusion coefficient, D, would be related to the temperature, T, using the Einstein relation where μ is defined as the mobility, k B is Boltzmann's constant, and T is the absolute temperature. By fixing T = 1, the diffusion coefficient would be directly related to the mobility. It is still interesting to notice that the mobility could, in principle, be directly interpreted as a sort of thermal parameter. In our model, the time period when the population is allowed to move without restrictions is characterized by a large value of σR (namely, σ high = 500 m), i.e., high temperature, while a lockdown is characterized by a lower σR (here, σ low = 0.5 × 10 0.5 m), i.e., low temperature. Thus, one can consider the infection rate problem in terms of heating/cooling of the system.

C. The infection probability
The core of the model, which contains most of the epidemiological data, is the probability for the ith healthy person to become sick. We assume that for each contact with another infected person, this process can be described by a Gaussian function of distance, weighted with a factor that parametrizes the sick person's conditions and social interaction, 17 where • ri(xi, yi) is the location of the ith healthy person and rj(xj, yj) is the location of the jth sick person, so |ri − rj| is the distance between them, • n sick is the total number of sick people in the area, • σr is the standard deviation (here, σr = 2.4 m) since recent studies show that even a slight breeze can drive droplets arising from a human cough over more than 6 m, 6 • f (a silent , nout) is a function that considers the social activity of the sick person and whether he has symptoms, which affect the spread of the virus outside the house. In our model, we estimate that infection by asymptomatic people is ∼50% lower than patients with symptoms, • ξ is a random number between 0 and 1, which allows us to consider some violation of the lockdown and the fact that even during a full lockdown, people continue to go out of their homes to buy grocery, go walking, etc.
In addition, we are assuming that each sick person will infect some of his household members. Since the latest estimates are that household infections are ∼15% from known cases (without lockdown 18 ), we estimated that the number of household infections is uniformly distributed between 0 and 3. This number is constant for all the simulations and does not depend on the population density or the lockdown constraint.

III. LOCKDOWN STRATEGIES
To date, lockdown has been imposed in many countries to reduce R 0 and, as a result, to increase T d , the doubling contagion time. It still remains unknown if the disease will start spreading again without control or if new local outbreaks will appear. In our work, we tested three different types of lockdown strategies. This choice is just representative of a potentially much broader set of options that could be analyzed using this method by simply varying the corresponding parameters and with a minimal computational cost. For all the models presented in this paper, we assumed the following conditions: • Days 1-14: no restrictions.
• Days 15-50: full lockdown with moderate social distancing (SD), i.e., people are forced to maintain a distance of 3 m from each other.
• Days 51-200: people must wear face-masks so that the daily infection probability (for the non-household members) is reduced to A very recent HKU hamster research study shows that by wearing a proper mask, the infection probability can be reduced by a factor of 3. 19 Therefore, given that not all of the population wears a mask properly and given the findings of Ref. 20, we estimated the probability of infection when wearing masks to be 1.4 times lower than that without masks.
Days 1-14 are the heating phase of the system. Thus, we expect the fastest increase in the number of patients these days. During days 15-50, we predict a phase transition from a hot system to a colder system (almost solid-like), which will reduce the infection rate. Each model was tested with the following constraints: • with and without moderate social distancing (SD) on days 51-200, • with and without infection from unknown sources (infected people not traceable to any known infection chain), • with and without strict SD for symptomatic patients after the 14th day (sick people with symptoms are forced to maintain a distance of 8 m from healthy people, equivalent to strict home isolation).
Hence, for each population density, N, we have 24 different simulations that we will use later to assess the effect of both SD and lockdown on the number of cases as a function of time.

IV. RESULTS
For the COVID-19 epidemic, the observed values of R 0 suggest that each infection directly generates 2-4 more infections in the absence of countermeasures like social distancing. 9,21 The doubling time, T d , is a function of R 0 22 such that the higher the R 0 , the lower the T d . In particular, T d = 2.5 days corresponds to R 0 = 4.
In our model, both R 0 and T d are directly obtained from the simulation and not pre-assumed. Figure 1 presents the doubling time, T d , as a function of the effective density, N. We have calculated the percentage of active cases from the total population in the first 14 days for each N. The calculation was performed in the heating phase only since there are no restrictions. This allows for a sort of calibration of the model to encompass the intrinsic features of the disease. Figure 1 shows that T d is a decreasing function of density. It was shown in Ref. 22 that for low R 0 (1 < R 0 < 1.5), the doubling time difference between high and low R 0 is reflected in our calculations by a high error for high doubling times, which are compatible with low population densities.

Physics of Fluids
The relation between T d and R 0 , as shown in Fig. 1, means that attempting to describe a wide area's current situation by means of some average value of R 0 might be highly inappropriate. On the other hand, the spread of the disease over smaller, more homogeneous areas in which the population shares a certain degree of mobility and social behavior would be quite well described by this parameter.
To give an example, a value of T d ≈ 2.5, corresponding to R 0 = 4, in the current model would correspond to N = 1.1 × 10 4 households per square kilometer. While this might appear as an unreasonably high density in an average urban context, it is still much lower than the average density in kindergartens, university classrooms, crowded social, religious, or sports events, etc. Hence, from now on, we will consider the value of N = 1.1 × 10 4 as representative of potentially dangerous situations present daily before the beginning of the pandemic. This will also show how different constraints affect the initial R 0 = 4, which is the highest estimation for R 0 without restrictions.
The various probability densities are sampled by means of standard techniques, in the spirit of a kinetic Monte Carlo simulation, for predicting the number of coronavirus cases as a function of time for different lockdown models and external constraints. Each case has been run 100 times. Note that each run starts with the same initial condition, only one sick person. Since the Monte Carlo algorithm is based on random numbers, we expect that every run will yield slightly different results. Repeating the simulation for 100 separated runs evaluates the algorithm's robustness. Results have been averaged and analyzed to determine the statistical error. Figure 2 shows our numerical results for the different simulations for the case of R 0 = 4 in the first 14 days. The upper (lower) panels of Fig. 2 are the numerical results without (with) unknown sources. In panels (b) and (d), the numerical results for the case where sick people with symptoms must maintain a distance of 8 m from healthy people (strict SD) are shown. As previously pointed out, the numerical results presented in Fig. 2 do not make any assumption on the doubling time, T d , but only on some observed features of the disease. Even though it is challenging to model each coronavirus infected area's specific characters, some characteristics are common to all simulations. The first 50 days have the same constraints-no restriction from the first day up to day 14 and a full lockdown from day 15 until day 50. Our numerical results show that even for R 0 = 4, forcing a lockdown after 14 days from the first case controls the spread of the coronavirus in a way that the peak number of cases occurs on day 24, while the rate of increase in the number of cases starts decreasing as of the 15th day. In addition, since the average time for recovery is 14 days (although up to 40 days for very severe cases) and the typical infection period ranges from the third day until the seventh day, the decreasing rate of the number of active cases during the lockdown is much slower than the increasing rate of the number of active cases without restrictions [as seen in various countries around the world (see Ref. 15 for up-todate data)]. Hence, we find that for an initial heating period of 14 days, the necessary lockdown period (i.e., the cooling time) required to reduce the number of active cases is much longer than 14 days.

Physics of Fluids
As of today, many countries examine different exit strategies due to the decreased number of active cases. In our simulations, we have tested several such exit strategies using different constraints. All of our numerical results indicate that the isolation of symptomatic patients [strict social distancing, panels (b) and (d) of Fig. 2] is effective and can reduce the peak number of active cases by a factor of about two without further restrictions, and up to a factor of 10 for model 3. In addition, from Fig. 2, we find that moderate social distancing can reduce the number of active cases but never as effectively as home isolation of symptomatic patients. Figure 3 is an inset of Fig. 2(c). From Fig. 3, is it easy to see the effect of the cyclic no-restriction/lockdown pattern. In principle, since the no-restriction/lockdown pattern is periodic in time, we would expect that the number of active cases as a function of time will have the same periodicity. For all the three models, there are no restrictions from day 51 until day 57 (green area, A). From day 58 until day 64 (red area, B), we impose a lockdown in models 2 and 3, which is reflected in a more moderate increase in the number of active cases than that of days 51 until 57.

ARTICLE scitation.org/journal/phf
is easy to see that the increased effectiveness of model 3 in reducing the number of active cases originates from the fact that, since it takes more than a week to cool the system (which is a result of the infection period), a week-week strategy cannot cause significant cooling of the system. An effective exit strategy might be based on time cycles and must include at least ten days of lockdown (see, for example, Ref. 9). In addition, one can use different (no time based) no-restoration/lockdown patterns such as those presented in Refs. 23 and 24. The lockdown strategies presented in Refs. 23 and 24 are not time-dependent; in contrast, they are based on optimization of a fixed number of infected/expired people over time. In general, our model can also be used for building these kinds of no-restriction/lockdown patterns for each population density.
In Fig. 4, we present our numerical results for a 4/10 day cyclic exit strategy for different σ low . Similar to Fig. 2, here, the upper (lower) panels are the numerical results without (with) unknown sources. For both cases, in the left panels [(b) and (d)], the numerical results are for the case where sick people with symptoms must maintain a distance of 8 m from healthy people (strict SD). Figure 4 shows that, as predicted, for instance, in Ref. 9, a 4/10 day cyclic exit strategy is useful for controlling the infection rate accompanied by home isolation of symptomatic patients. The comparison between Figs. 2 and 4 indicates that although the maximal percentage of active cases is similar for both the week/two-week pattern and the 4/10 pattern, there are differences between the two patterns. For the 4/10 pattern, this cyclic pattern induces a moderate increase in the percentage of active cases but without a local decrease in the number of patients until the peak on the 130th day. In contrast, for the case of the week/two-week pattern, the two-week lockdown will cause a local decrease in the percentage of active cases and a more steep increase in the percentage of active cases during the no-restriction week accomplished by a peak in the percentage of active cases on the 110th day. Hence, the comparison between the two patterns implies that a cyclic no-restriction/lockdown pattern can control the spread of the epidemic in daily life, even for an initial doubling time of 2.5 days.

V. COMPARISON WITH REAL DATA (SWEDEN)
From the beginning of March 2020, Sweden took a different approach from the rest of the world by not imposing a policy of lockdown on its citizens. Therefore, it is of interest to examine the rate of increase in the total number of active cases in Sweden from the beginning of March until today (see Ref. 15 for up-to-date data), compared to our model under the assumption of no restrictions. In Fig. 5, we show the total number of active cases in Sweden (normalized) from the 100th case until today, in comparison to our predictions for a population density of N = 3500 households. Note that this population density is much more diluted than the population density used for the previous simulations, which may explain the slow rate of increase in the virus spread even when there are no limits. Figure 5 shows that our model can predict the spread of the virus for different societies, reflected with varying densities of population.

VI. CONCLUSIONS
This paper presented a kinetic Monte Carlo algorithm for modeling different scenarios of the infection rate of the novel coronavirus disease. This model's main feature lies in its extreme flexibility and in the fact that the parameter R 0 is obtained from the simulation and not pre-assumed. It can rather be used, in principle, as a way to tune up the other parameters better based on the post-processing of clinical and epidemiological data.
Although it is challenging to model the specific characters of each coronavirus infected area, our results show that strict social distancing and a cyclic time pattern might help to keep the infection rate under control over a long period, even for an intrinsic doubling time of 2.5 days and in the presence of infection from unknown sources. Our ability to model and prove differences between the different lockdown patterns sharpens the need for physical and mathematical models that allow examining different ways for reducing the spread of the epidemic. From the physical point of view, effective strategies for controlling the infection rate of a specific area should lower its effective temperature as much as possible by keeping social distancing and avoiding creating hot spots such as those related to high concentrations of people on a daily basis.