Monte Carlo simulation of COVID-19 pandemic using Planck's probability distribution.

We present a Monte Carlo simulation model of an epidemic spread inspired on physics variables such as temperature, cross section and interaction range, which considers the Plank distribution of photons in the black body radiation to describe the mobility of individuals. The model consists of a lattice of cells that can be in four different states: susceptible, infected, recovered or death. An infected cell can transmit the disease to any other susceptible cell within some random range R. The transmission mechanism follows the physics laws for the interaction between a particle and a target. Each infected particle affects the interaction region a number n of times, according to its energy. The number of interactions is proportional to the interaction cross section σ and to the target surface density ρ. The discrete energy follows a Planck distribution law, which depends on the temperature T of the system. For any interaction, infection, recovery and death probabilities are applied. We investigate the results of viral transmission for different sets of parameters and compare them with available COVID-19 data. The parameters of the model can be made time dependent in order to consider, for instance, the effects of lockdown in the middle of the pandemic.

Introduction

The COVID-19 pandemic, that has so far caused over six million deaths worldwide (Hui et al., 2020, Anon, 2022a), has recently induced numerous studies on the mathematical modeling of temporal distributions of infected cases and fatalities (Garin et al., 2020, Cooper et al., 2020b, Cooper et al., 2020a, Langel, 2021, Purkayastha et al., 2021). Various statistical models for pandemic spread and their predictive power have been benchmarked against data (Flanders and Kleinbaum, 1995, Amaro et al., 2021).

The SIR (susceptible–infected–recovered) model – developed by Kermack and McKendrick (1927) – and derived models are the most widely used to study viral spreading of contagious epidemics or mass immunization planning (Weiss, 2013, Chauhan et al., 2014, Chao and Dimitrov, 2016) since they provide mean values of the cumulative incidence as a function of time in a deterministic way. SIR models belong to the compartmental type, where the equations only involve the time variable, and the incidence functions refer to the total number of individuals in each compartment or subset of the total number of individuals in the system.

Generalization to more realistic models of the infection transmission process requires solving spatio-temporal equations (Ganesan and Subramani, 2021) in a two-dimensional space representing the surface of a region under epidemics. Such individual-level models were first studied by Kendall (1957) and  Bartlett, 1957, Bartlett, 1956, who also considered both deterministic and stochastic kind of descriptions (Allen, 2008, Andersson and Britton, 2000a).

Some stochastic approaches are based on SIR epidemic models with stochastic perturbations (Rao, 2014). The SIS model allows the infectivity to fluctuate around a mean value by introducing external fluctuations in the SIS epidemic model (Otunuga, 2019). Another example is the stochastic SIRC epidemic model with discrete time delay (Rihan et al., 2020).

Stochastic epidemic models can also be formulated in terms of continuous-time Markov chains (CTMC), stochastic differential equations (Allen, 2017), multivariate stochastic processes (Kim et al., 2018), or Brownian motion (Ming et al., 2016). A recent CTMC study of the transmission dynamics of the COVID-19 in Wuhan, China was presented in Olabode et al. (2021) (see Ref. Andersson and Britton (2000b) and references therein for a general description of stochastic epidemic models and methods for their statistical analysis).

In this paper we present a study of space–time propagation of a viral infection using a new stochastic model inspired by concepts of statistical physics and nuclear physics. This model is compared to a simple compartment model developed in Ref. Amaro et al. (2021) (the D model), where we studied the accumulated fatalities for a series of countries during the first COVID-19 wave.

The purpose of this work is to present the Monte Carlo Planck (MCP) model in order to study the stochastic behavior of the spatio-temporal propagation, to test the predictive power of these kind of models with respect to COVID-19 available data as well as to assess the validity of the D model.

Further Monte Carlo-type studies have been developed to simulate the pandemic transmission (Xie, 2020, Maltezos and Georgakopoulou, 2021), but use other approaches such as, for instance, a 2D-random walk Monte Carlo simulation based on proximity infection spread. In the present MCP model we take a physical point of view, where the infection is a result of an interaction, and it follows to some extent statistical laws similar to the interaction in a system of many particles. Thus we model the infection as a likely result of a physical interaction. We apply concepts taken from statistical and nuclear physics such as interaction range, energy, temperature, and interaction cross section. These concepts, adapted to the present case of interest, are useful to describe probability distributions and interaction coefficients in physical systems, which are here characterized by epidemiological probabilities of infection, recovery and death.

Therefore, the novelty of the present approach with respect to previous stochastic models is that we describe the spread of the disease using magnitudes and concepts taken from physics. Our approach is that the epidemic occurs in two steps: i) by interaction between individuals in a population, followed by (ii) a probability of infection subsequent to the interaction. Without interaction there would be no infection.

In our model the population would be similar to a low-density gas with molecules that interact through elastic or inelastic scattering that can change the internal state of the molecule (infection). A fundamental microscopic variable that determines the epidemic is the number of interactions for each particle per day ( day is the chosen time step). This number is computed randomly for each particle of the system and for each time interval . In our model, the value of is determined by three parameters, temperature, range and cross section.

In physics, the temperature of a system determines the average value of the energy of a particle. Particles with more kinetic energy move faster and are more likely to interact. Thus in the case of a population, we define energy as a number, , describing the mobility of individuals. People who often go out () for a walk, shopping or for work will be more likely to interact than people who rarely leave home (). In our model the energy, , follows the same probability distribution as the number of photons in a radiant cavity with temperature (Planck’s law) as described in the next section.

The range is a geometrical parameter giving the average distance that the particles move away from their equilibrium position. People who move within a larger range can spread the infection over greater distances, for example if they travel often, as opposed to people who stay in the vicinity of the home.

The cross section in physics is proportional to the probability of interaction between two particles. For example, when a photon strikes an atom or molecule there is a probability of interaction by absorption or scattering, or it may pass through without interaction. The concept of cross section can be illustrated with the following example. Imagine an experiment in which we randomly throw darts over a target of cross section . If the beam of projectiles spans a cross-sectional area . The probability that a dart hits the target (interaction probability) is and the expected number of darts hitting the target is Similarly, in the case of a population we will characterize the probability of interaction between individuals by a cross section . The interaction probability quantifies in our model the exposure to the disease, which depends, for example, on protective measures such as keeping a distance, wearing gloves and a mask, etc.

The three “physical” parameters of the model, , and , together with the “epidemiological” parameters, probability of infection, recovery and death, will be fitted to the data of daily deaths in several countries.

The paper is organized as follows: Section 2 introduces the MCP model. In Section 3 we present results of the simulation applied to COVID-19 fatalities data for several countries. Finally in Section 4 we draw our conclusions.

The Monte Carlo Planck model

We consider the system of a two-dimensional lattice or grid, where each unit cell has two coordinates , with , where is the total size of the system, with surface . Each cell is occupied by an individual, who can be in any of four states or classes: susceptible, infected, recovered or dead. This state is specified by four fields or matrices , , , and that can take the values 1 or 0 (state quantum numbers), depending whether the individual in cell belongs or not to the corresponding class. The model can be easily extended to include other classes or intermediate states.

In the initial state, for , all the individuals are susceptible. Then for all except for one initially infected cell , where We choose the initial infected cell at the center of the grid, but in the model it can be chosen randomly or even to be several infected cells.

In the simulation we compute the state of the system in time intervals of day. At the end of the day, , we compute and store in the matrices , , , the current state of the system. We also store the total number of susceptible , infected, , recovered, and dead, and the daily increments, , , , , and repeat the calculation for the next day.

In the simulation for day we apply the following algorithm to each cell .

If it is not infected, , do nothing.

It it is infected, , then we assume that it can infect only people in other cells within a finite region of size . The value of the range takes into account the zone of influence of individuals in one day. To simplify the algorithm in Cartesian coordinates we assume that the region is a square of side , centered at the cell , but it could be also circular or other shape. Now the Monte Carlo starts by computing the value of randomly with some probability distribution. We assume an exponential law which indicates that is less probable to move far away from the average interaction range of the individuals. This is one of the parameters of the model.

Once is chosen for an infected particle, it can interact randomly with any of the individuals within the region with area . As in particle physics, we characterize the interaction between individuals by a cross section . The number of possible interactions is then given by the probability formula when particles are shot over a target containing particles (see Eq. (1)), where is the surface density and the number of shots is related with the mobility of the individual within its zone of influence, indicating if it has a narrow social behavior and how often it moves around each day. We call this property the energy of the individuals. For instance, Young people have more energy than old people, or people with a job where they interact often with other people. We compute the value of the energy randomly with an energy distribution. In analogy to the Planck law for photons from statistical physics, we write the energy distribution in a simplified form as where the parameter is the temperature of the system that determines the average energy or number of shots. Eq. (6) arises from Planck’s hypothesis for the probability of an electromagnetic wave with frequency to have quantized energies of . where is determined by demanding that In our model Plank’s constant, , Boltzman’s constant, , and the frequency, , are not relevant and we have taken them equal to one. Thus for each infected cell its energy is obtained as a random number generated with the probability distribution in Eq. (6). In this work we choose unit cells and therefore assume in arbitrary units, but the model can be modified to include specific densities and sample population sizes with known number of individuals. In the model the information on the density is included in the range . To increase the density is equivalent to increase the range and the temperature. The result of this work are given for unit density and the values of and should be appropriately re-scaled to apply the results to different densities.

Once we know the range and the number of interactions of an infected cell, we choose randomly cells, with coordinates , within the interaction region. If cell is susceptible, , it can be infected with probability . We then decide randomly it it is infected or not according to that probability. In case that it becomes infected, we change the value of the corresponding matrix elements and . We also store the time it becomes infected in the infection time matrix . This matrix is important for the next step in deciding when the particle becomes recovered or dead.

After interacting with all the individuals within the range , along the day, finally we decide if cell becomes recovered or dead at night. The removal probability is the probability of recovery or death. It must be a function of time with the condition that it is very small for and is equal to 1 for large t. In our model we use the following function that verify these conditions: where is the time such that and , so that the removal probability is very small for low times, and later it increases with time at a rate driven by the parameter . For , the removal probability is one. All the individuals recover or die. Therefore we randomly decide if is removed at the end of the day with a probability that depend on the time that particle has been infected. In case is removed, we change the state quantum number in the infected matrix, .

In the final step, after particle is removed we randomly decide whether it dies or recover. The probability of death is a parameter with . In the death case the function is increased by one. Otherwise it is recovered. Finally we store the matrices (0), and (1) accordingly.

We repeat the procedure on the next cell.

The algorithm so defined depends on the following parameters

The size of the system , The unit of length in natural units, is the average distance between individuals. This determines the number of cells .

The range of interaction , indicating the maximum distance on average that the infection can propagate in one day. This parameter can be time dependent, and, by changing it, one can simulate for instance lockdown or long travels in holidays.

The cross section measures the probability of interaction between two individuals. In the case of classical particles interacting with contact forces the cross sections is just the effective geometrical area of the pair. For long range forces the cross section is larger. In our case it is typically much larger than the geometrical human size, because this parameter takes into account the human behavior. Humans move and are social increasing its cross section, like a long range force.

The temperature of the system, determines the average energy of the individuals, which is related to the mobility within its range of interaction.

The infection probability

The removal probability parameters and (two parameters).

The death probability of removed

These are in total 8 parameters , , , , , , and . In addition, each one of these parameters can be easily made time dependent to study the effects of political enforcement, lockdown, social distancing, large events, etc. The model could also be modified easily to include coordinate dependence of the parameters to specify for instance a large event in some region by increasing the temperature or the cross section in that region.

Actually one can argue that increasing temperature, increasing cross section and increasing the infection probability are expected to produce similar effects on the pandemic evolution, and that the pandemic evolution is effectively driven by less parameters, perhaps six or less, unless time or space dependence (or non constant densities) is introduced in them.

In the next section we fit the parameters of this MCP model to describe data from the first wave of the COVID-19 pandemic.

Results

Firstly, the MCP model requires fixing the number of individuals, which is equal to the area of the square lattice (each individual occupies one square). To make the calculation manageable, we chose an initial size of , that is individuals in the sample, although we latter also consider the effect of making L larger. Obviously, if the number of deaths exceeds this figure, to study a more realistic case, a size comparable to the population of the country in question should be considered. But in this first simulation we prefer to make the model as simple and manageable as possible in order to investigate whether under these simplifications it is still possible to describe the data.

An important detail is that, given the probability of death , the expected mean number of deaths is fixed from the beginning . If this value is not known experimentally, we first fix it by fitting a SIR-type model to the incomplete data. Specifically, we apply the D- and D2-models of Ref. Amaro et al. (2021), which consist of sums of logistic functions, to estimate the probability of death before doing the Monte Carlo. The death number as a function of time in the D-model is where the three parameters are fitted to data and the initial time is arbitrary. In the case of the D2-model we fit the sum of two D-functions with six parameters in total. This provides a method for estimating a prediction for the total number of deaths ( in the case of the D2-model). This value is used as input to set the parameters of the Monte Carlo simulation since we know that the D model works well as a starting point for pandemic forecasting (Amaro et al., 2021). Thus the Monte Carlo will be useful to study the stochastic distribution of daily cases, since the SIR models only provide the mean values without statistical fluctuations.

The simulation was started assuming that in the initial instant there is an infected individual in the central cell of the array and all the others cells are susceptible. The algorithm of the previous sections is then applied by assuming time steps day. In Fig. 1 we show the results of the Monte Carlo simulation of the deaths-per-day, , in Spain. We fitted the data until May 30, 2020, where clearly the first COVID-19 wave had practically ended. We first fitted the D2 model (as the sum of two logistic functions) to the data to better estimate the end of the epidemic, and then we fitted the MCP parameters to the curve of total accumulated deaths of the D2-model. It is worth mentioning that fitting the Monte Carlo to the accumulated deaths curve is more reliable than fitting the daily data, , and gives better results. This is because the cumulative curve data has less fluctuations and less error than the daily data.

Fig. 1

Fit of the P-model to the cumulative deaths and deaths per day by COVID-19 in Spain during the first wave. The red points are the data compiled up to May 25, 2020. The two yellow lines in the left panel are the two D-functions appearing in the D2 model, Eq. (10).

The MCP parameters are given in Table 1. The parameters of Table 1 have been tuned to describe the experimental data, including the lockdown. In general The Monte Carlo points have to be shifted in time ( days in the Spain case) to agree with the center of the experimental peak, that has a different time origin. This is so because there is an arbitrariness regarding the start of the pandemic.

Table 1

Parameters of the MCP model for Spain, Sweden and South Africa.

	L	R_0	ρσ	T	t_r(days)	b(days)	p_i	p_d
Spain
First wave	300	8 (3)	0.1	20	20.7	4.5	0.2	0.3
Sweden
First wave	300	10	0.3	20 (15)	24	4	0.2	0.3
South Africa
First wave	300	2.9	0.09	20	20	4	0.15	0.2
Second wave	300	3	0.1	20	20	4	0.15	0.5

The value of the range during lockdown in Spain and the temperature during gathering restrictions in Sweden are in parentheses.

For Spain, the MCP simulation includes lockdown effect on day 20, where the value of the range was reduced to (fitted also to the data). In Spain the lockdown started on day 14M (March 14), corresponding to day 7 in the plot of Fig. 1, and the maximum of the daily deaths, , occurred on 2 A (April 2), corresponding to day 26 in the plot. The distance between maximum and lockdown therefore is days. In the MCP simulation this distance is 15 days.

From the results of Fig. 1 we see that the MCP and the D2 models present the same trend. Both models describe the death data, both total and daily deaths, very similarly. This means that the hypothesis of the MCP model that the spread of infection can be described by physical interactions in a system of many particles is correct, since it agrees with the statistical model of the SIR type and is able to describe the experimental data of the pandemic.

The main difference with the D2 statistical model is that the MCP acquires random statistical fluctuations, which the SIR models do not have, since they only describe mean values. We also see that the fluctuations in the data are much larger than in our model, even considering the statistical errors of the daily data. These values have been taken from the Poisson distribution . However it must be taken into account that the official death data have systematic errors due to the data collection methods and delays, and that these errors could be quite large. In such cases it may be reasonable to first compute a moving average of the data over three or more days to reduce fluctuations before fitting the models.

Concerning the values of the parameters, the death probability is fixed by the total fatalities. In Spain the three physical parameters are the following. The initial range , later changed to after the first week to take into account the lockdown. The temperature or mobility , means that in average 20 shots are used to compute the interaction in Eq. (4). The cross section multiplied by the density gives the probability of interaction between cells . This means about 2 interactions in average. The infection probability is . Finally the recovery time is days and the recovery interval days These parameters give an idea of how long the infection lasts, until there is recovery or death, which is around three weeks.

In Fig. 2 we show a simulation of the case of Sweden. In order to keep the probability of death equal to the Spanish one, we have randomly chosen 85% of the cells with total lockdown, which produces the same effect as reducing the size of the susceptible ones, in order to obtain the number of deaths from the first or . As a result, the value obtained for the range is larger than in Spain , and also the value of . In Sweden there was no lockdown and this value does not change in time. But in Sweden there was a gathering restriction that no more than 10 people could meet. This has been simulated in the Monte Carlo by lowering the temperature from to during gathering restrictions. The gathering day here has been taken equal to 40. The recovery time is about three days larger than in Spain.

Fig. 2

Fit of the P-model to the cumulative deaths and deaths per day by COVID-19 in Sweden during the first wave.

In Fig. 3 we present results in the case of South Africa. The first and second waves are well differentiated and each of them can be described separately with model D. This is what has been done in the figure, adjusting the data for each wave separately. In the same way, the MCP model has been adjusted here separately for each of the two waves. The first Monte Carlo was adjusted for data through August 30, 2020. The second Monte Carlo was adjusted for data from the second from this date through February 10, 2021. Each Monte Carlo simulation has been adjusted independently with a lattice of . That is why the probability of death in the first wave is 0.2 whereas in the second wave is 0.5, since the number of deaths according to these models would be in a 5:2 ratio between the second and first waves. The remaining parameters of the Monte Carlo have not changed from those of the first wave, except for one decimal in the case of the range and the exposition . The right panel in Fig. 3 is a Monte Carlo prediction of the second wave daily deaths from the data up to February 10, 2021.

Fig. 3

Fit of MCP model to the cumulative deaths by coronavirus in South Africa during the first (left) and second (middle) pandemic waves, together with the overall daily deaths (right). Data from (Anon, 2022a) fitted up to Feb. 10 2021.

Note the similarity between the parameters of the model for the three countries that are shown in Table 1. This shows some universality in the physical characteristics that describe the process in our simplistic model. A more realistic model would require describing the geographic details and the average population distribution according to a more realistic geometry than that taken on a square lattice with constant density. But in view of the results obtained there is hope that this model will be useful to describe a more realistic case, but requiring some more work.

To finish, we have also developed (Anon, 2021c) an educative android app to demonstrate the virus propagation interactively by plotting the daily grid evolution (Fig. 4) using a color for each cell state: susceptible (green), infected (red), recovered (blue) and dead (black). The resulting SIR functions S(t), I(t) and R(t) computed with the MCP model are also displayed in the app screens, as well as the daily deaths and cumulative deaths. The initial pandemic parameters are given by the user in the initial screen of the app. The Monte Carlo is started by pressing the running button and the grid evolution is shown day by day on the screen. Here the user is allowed to change interactively the parameter range, , exposition, , and mobility, , to observe the effect of contention measures and stop the pandemic from spreading. More information, including online lectures and the Fortran code, can be found at UWC’s nuclear GitHub (Anon, 2022d).

Fig. 4

Monte Carlo Plank model simulator for Androids illustrating how the pandemic spreads in a lattice of length – which corresponds to a population of – as a function of time.

Conclusions

In this work we have presented a Monte Carlo model (MCP) of the spread of an infection, based on interaction parameters inspired by the physics of many particles. Starting from physical quantities such as interaction range, cross section and temperature, we translate these concepts to the mathematical description of a pandemic under the names of range, exposition and mobility.

The model consists of individuals in a rectangular lattice that can interact with their neighbors within a range and at a certain speed, with a certain probability of interaction, which in turn can produce infection with the result of recovery or death.

The model fits only for daily deaths, although it could be adapted to fit other infection data. In this work we have proceeded to calculate the simplest case with certain simplifications to keep the calculation as simple and short as possible. In particular, we have set a square lattice with 90,000 individuals, with which the probability of death is fixed by the data when they are adjusted to a SIR-type model.

In general, the model is capable of describing the experimental data as well as the SIR-type models, with the difference that the MCP has a stochastic behavior, as it is a Monte Carlo simulation subject to probabilistic parameters. It also allows the inclusion of parameter modifications over time to simulate mobility restrictions or quarantines.

The model parameters have been adjusted to describe pandemic waves in Spain, Sweden and South Africa, with similar results for the parameters, showing a certain universality of the physical process involved in the pandemic.

In future work we will study the dependence of the results with the parameters of the model, as well as extensions to describe larger populations, of the order of magnitude of complete countries to determine parameters in more realistic cases.

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.