A random walk Monte Carlo simulation study of COVID-19-like infection spread
aa r X i v : . [ q - b i o . P E ] J un A random walk Monte Carlo simulation study of COVID-19-like infection spread
S. Triambak ∗ and D. P. Mahapatra † Department of Physics and Astronomy, University of the Western Cape, P/B X17, Bellville 7535, South Africa Department of Physics, Utkal University, Vani Vihar, Bhubaneshwar 751004, India (Dated: June 30, 2020)Recent analysis of COVID-19 data from China showed that the number of confirmed cases followeda subexponential power-law increase, with a growth exponent of around 2.2. The power-law behaviorwas attributed to a combination of effective containment and mitigation measures employed aswell as behavioral changes by the population. In this work, we report a random walk MonteCarlo simulation study of proximity-based infection spread. Social distancing is incorporated in thesimulations through a single parameter, the size of each step in the random walk process. The stepsize l is taken to be a multiple of h r i , which is the average separation between individuals. Threetemporal growth regimes (quadratic, intermediate power-law and exponential) are shown to emergenaturally from our simulations. For l = h r i , we get intermediate power-law growth exponents thatare in general agreement with available data from China. On the other hand, we obtain a quadraticgrowth for smaller step sizes l . h r i /
2, while for large l the growth is found to be exponential.Together with available data, these results suggest that the early containment of the disease withinChina was close to optimal. We further performed a comparative case study of data from threeother countries, India, Brazil and South Africa. We show that reasonable agreement with thesedata can be obtained by incorporating small-world-like connections in our simulations. Introduction
Following its outbreak in the Hubeiprovince of China, the global spread of the novel coron-avirus disease (COVID-19) has reignited efforts to betterunderstand infection spread and mortality rates duringthe pandemic. Significant emphasis was placed on mod-eling the spatio-temporal spread [1–3] of the disease, inorder to make reliable predictions. A key statistic in suchepidemiological analysis is the basic reproduction num-ber R , which defines the expected number of secondarycases from one infected individual in a completely sus-ceptible population. Data from the early phase of theCOVID-19 outbreak showed good agreement with mod-els that assumed an exponential growth of infections intime ( t ), with a mean R ranging from 2 .
24 to 3 .
58 [4, 5].However, subsequent laboratory confirmed cases in Hubeishowed that soon after the initial stage, the temporalgrowth in the cumulative number of infections ( N ) wasinstead subexponential and agreed reasonably well witha power-law scaling N ∝ t α [6]. This was consistent withdata from other affected regions in mainland China (with α = 2 . ± .
3) and was attributed to a depletion of thesusceptible population due to effective containment andmitigation strategies that were put in place and followedafter the initial unhindered outbreak [6].A potential stumbling block in such analyses is thatthe reported number of infected cases may be inaccu-rate, due to a non-uniform sampling of the entire sus-ceptible population in a given region. In such a scenarioone can alternatively examine the number of reporteddeaths (due to COVID-19 complications) as a function ∗ [email protected] † [email protected] of time. This is justified, as the number of deaths aregenerally more accurately recorded and are always a frac-tion of the total infected population. Indeed, mortalitydata from the National Health Commission of the Peo-ple’s Republic of China and Health Commission of theHubei Province showed similar power-law behavior, withan exponent α ≈ . Monte Carlo simulations
Our simulations are basedon minimal prior assumptions and use a two-dimensionalrandom walk model, within which the individuals from N i n f e c t ed α = 2.24 ± 0.05 α = 2.27 ± 0.06 α = 2.24 ± 0.08 Number of time steps (t)Population10k Population 6k Population 2.5k
FIG. 1. Simulation results for three population sizes with a fixed number density of 10k/unit area. The simulations assumethat all random walkers take steps of length h r i , their mean separation distance. The dashed lines are best-fit results for power-law growth N = At α . The quoted uncertainty for each extracted α includes a ± σ statistical uncertainty and a systematiccontribution. The latter were estimated by applying a conservative one channel shift to the data along the time axis andredoing the fits. an entire susceptible population are independent randomwalkers. They are initially represented by uniformly dis-tributed points in a given isolated region. The simu-lations begin with an initial condition of one infectedwalker, assuming all other points are ‘normal’, yet 100%susceptible. As the simulation progresses, each walkertakes steps in random directions, defined by a suitablestep size described below. The ‘spread’ of the disease issaid to occur whenever an infected point comes within a‘touching’ distance from a ‘normal’ point, thereby pass-ing on the disease. This assumption includes indirecttransmissions through infected surfaces, etc . An impor-tant aspect of the simulations is that in addition to theirspatial mobility, the incremental number of simultaneousdiscrete steps taken by all random walkers quantifies timeprogression. This allowed us to track the evolution of theinfection spread, in terms of growth in the number of in-fected points with respect to the number of ‘time-steps’( t ). Our simulations also incorporate the possibility of having onefraction of the population to immune to infection, while allowinganother fraction to recover from the disease, with a possibility ofreinfection. However we do not delve into these aspects here asthey are beyond the scope of the present work.
To put the above in perspective, for N randomly dis-tributed points over area A , the mean distance of sepa-ration h r i between any two random walkers is ∼ p A/N .Therefore, for a metropolis such as New York city, whichhas a population density of ∼ h r i is ∼
10 m. For an arbitrarily sized region, given the fixedpopulation density, this would result in 10,000 points perunit area, with h r i = 10 − length units and a maximum‘touching’ separation of 2 × − units. The latter corre-sponds to ≤ , assuming that eachwalker’s mobility is constrained. This was imposed bythe condition that all members of the populace only takesteps of length l = h r i , the average distance betweenthem. Effectively, l quantifies the social distancing aspectin the simulation. The l = h r i condition ensured that All simulations described here were carried out for N randomwalkers over a unit area.
200 400 600 800 1000
Number of time steps (t) N i n f e c t ed l = 5 〈 r 〉 l = 〈 r 〉 l = 〈 r 〉/2 l = 〈 r 〉 /3l = 〈 r 〉 /4 N = . e . t N = . t . N = . t . N = . t . N = . t . FIG. 2. Simulated growth in the number of infections andtheir best-fit curves, obtained for different step sizes taken bythe random walkers in a population of 2.5k/unit area. Be-low a threshold value of l ∼ h r i /2, the growth is observed tobe nearly quadratic, regardless of the step size. The otherextreme shows exponential behavior, while power-law growthsimilar to what was observed in China [6] lies in the interme-diate regime. the infection could only spread through ‘nearest neigh-bor’ interactions and that each random walker is on aver-age confined within a local neighborhood. We performedthree such simulations for a fixed density of 10k/unit areaand three population sizes 10k, 6k and 2.5k respectively.Simulation set II (Fixed population and density, differ-ent step lengths): In the next step we probed the depen-dence on both population density and l in five separatesimulations. These simulations assumed a smaller den-sity of 2500 walkers/unit area, and different step sizesfor the walkers, with lengths h r i / h r i / h r i / h r i and5 h r i . The results are discussed below and in Fig. 2. Results and analysis
In Fig. 1 we plot the growthin the cumulative number of infected points, obtainedfrom simulation set I. The results show that independentof population size, the number of infections follow a t α power-law growth in time, with α about 2.2. Whilethe power-law behavior may not be completely unex-pected [13], it is interesting that we obtain very similarvalues of near-quadratic exponents, as observed with thedata reported in Refs. [6, 7]. In simulation set II, forstep length h r i , we determine almost identical power-lawgrowth as in Fig. 1, again in agreement with the observa-tions of Refs. [6, 7]. This is shown in Fig. 2. Addition-ally, we also find that our extracted power-law exponentsare consistently similar for this step size, regardless ofthe population density used. In comparison, if all mem-bers of the sample population were to take larger random We do not quote uncertainties in the fit parameters here, as thisfigure only serves to highlight the systematic trends of the curves. steps of length 5 h r i , on average interacting with pointslocated much further away than their nearest neighbors,we find that the number of infected individuals blows uprapidly, showing near exponential behavior. This wouldbe similar to a scenario where no control interventions arein place or being followed. Not surprisingly, unlike themore constrained random walk, the slope for exponentialgrowth is found to strongly depend on the populationdensity, and is larger at higher densities. Figure 2 alsoshows the other extreme in terms of the temporal growth,obtained using step sizes smaller than h r i . As apparentin the figure, the results from these simulations showedquadratic growth, for all step sizes less than a thresholdvalue of around h r i /2. This effectively implies a lower-bound on the growth exponent ( α = 2), exactly as in thecase of peripheral spreading [8]. Furthermore, the effectof social-distancing is clearly evident from the observeddelay in reaching the saturation value. A few importanttakeaway points from the above analysis are: i) A step-length parameter of size h r i reproduces the growth curvesfrom China reported in Refs. [6, 7]. ii) Under such con-dition, the growth exponent is found to be independentof the population density. iii) On further restricting therandom walkers’ mobility by reducing their step sizes, weobserve that the slowest unavoidable temporal growth isquadratic in nature. This is similar to peripheral growth,proposed in Ref. [8]. iv) Since quadratic scaling appearsto be a limiting case, it seems unrealistic for a large pop-ulation to achieve such growth. Therefore, at face valuethe available data suggest that the containment measuresand response in the eight affected Chinese provinces men-tioned in [6] were most likely optimal and could not havebeen significantly improved upon any further. Power-law, exponential growth and small-world-likeconnections in observed data (India, Brazil and SouthAfrica)
To make further comparisons, we examine thedaily growth in the cumulative number of deaths (whichis a fraction of the total infected population) reported [12]for three countries, South Africa, India and Brazil. TheseBRICS countries were not arbitrarily chosen, as theypresent an interesting comparative study for several rea-sons. They face common challenges in terms of povertyand economic inequality, and are home to some of themost densely populated informal settlements in the world(such as Khayelitsha in Cape Town, Dharavi in Mumbaiand the favelas of Rio de Janeiro and S˜ao Paulo). Thelack of proper sanitation is a common theme in these im-poverished city pockets, where, given the circumstances,expecting the residents to follow strict social-distancingprotocols is a tall order [14–16]. Secondly, the responseof the political leadership of Brazil to the COVID-19 cri-sis was strikingly different from the governments of Indiaand South Africa. While the latter two countries swiftlyimposed extended periods of severe lockdown [17–20] be-ginning in the month of March, Brazil has thus far notpursued a concerted policy for such containment [21]. As
10 20 30 40 50 60 70 80
Days starting March 20 C u m u l a t i v e nu m be r o f dea t h s
10 20 30 40 50 60 70
Days starting March 28
Days starting March 18
10 50 0.1k1k2k5k 10 50 0.1k1k10k10 30 50 70200400600
India BrazilSouth Africa
N = A t α with α = α = α with Q u a d r a t i c E x ponen t i a l N = A e α x with α = FIG. 3. Cumulative number of deaths reported for India, South Africa and Brazil from the WHO, until June 1 [12]. Thefitted growth curves are shown together with their 95% C.L. bands. The Brazil and India data show a single growth-exponent,consistent with power-law behavior. On the other hand the overall data for South Africa suggests exponential growth. Closerinspection shows that the initial growth phase for South Africa was quadratic. The growth exponents are quoted in the insetswith ± σ uncertainty. of today Brazil is poised to become the epicenter of theCOVID-19 pandemic. The cumulative death data re-ported for the three countries, with their correspondingfits are plotted in Fig. 3. While we do observe power-lawgrowths with exponents of of 2 . ± .
02 and 3 . ± . Number of time steps (t) N i n f e c t ed
3% of the population have step size 5 〈 r 〉 No small world connections
São Paulo (7.2k/Sq. km.) α = . ± . α = . ± . FIG. 4. Comparison of power-law growth exponents for anexample city such as S˜ao Paulo (population density 7.2k/unitarea), with and without small-world-like connections in thesimulation. The latter were introduced by having a random3% of the population take steps of length 5 h r i . The uncer-tainties in α are quoted similarly as in Fig. 1. in the population (representing essential service workersand non-compliant citizens etc.), who are allowed to takemuch larger randomized steps, bounded by the area A .If strict containment measures were not adhered toby members of the population, it would correspond toa combination of two effects in our random walk model:(i) All the random walkers would travel longer distancesthrough larger step sizes. (ii) A small fraction of the pop-ulation has much longer-ranged mobility compared to theabove. This establishes small world connections betweeninfected individuals and the rest of the susceptible pop-ulation.Our simulation results for an example city such as S˜aoPaulo (with a population density of 7,200/unit area) areshown in Fig. 4 For a uniform step length h r i we deter-mine a power-law exponent of 2 . ± .
06, again in agree-ment with our previous observations. On increasing thestep lengths of a randomized ensemble comprising 3% ofthe city’s population to l = 5 h r i , we find that the expo-nent increases to 3 . ± .
08, very similar to the Brazildata shown in Fig. 3. This clearly shows the contributionof long distance movers to the spread of the pandemic. Asa matter of fact, it is well known that long-ranged travelcan dramatically accelerate the spread of infection [24].The higher exponent for the data from India can beexplained similarly. Despite its best attempts, the coun-try’s COVID-19 containment strategy was challenged bythe sheer scale and diversity of its population. For exam-ple, it is known that on several occasions people defiedsocial-distancing measures to attend religious gatheringsin large numbers [25, 26]. Furthermore, the sudden andunprecedented lockdown in India resulted in a humani-tarian crisis, with millions of daily wage inter-state mi-grant workers from the rural hinterland left jobless inthe big cities [27]. This led to a large-scale migrationback home for thousands of such families, many of themtraversing large distances of the country on foot [28, 29].The exodus has only intensified with the relaxation ofrestrictions in May. It is widely expected that the cu-mulative growth of infections in India will continue toincrease as the lockdown ended around June 1 [30]. Asimilar rise is also expected for South Africa, as it movedinto a much more relaxed lockdown scenario (Level 3) onJune 1.Recently, there have been several attempts tofit the sigmoid-type curves for country specificCOVID-19 infections, with a generalized logistic func-tion (see Ref. [31] and references therein) of type N ( t ) = A (1 + Be − µt ) / (1 − m ) , that solves the Richardsdifferential equation [32]. We caution that such an ap-proach can lead to inaccuracies, particularly when a con-tainment policy is followed. It is clear that the expres-sion for N ( t ) does not produce a linear relationship be-tween ln N and ln t , as expected for (contained) power-law growth. This is manifest in the results of our simu-lations and validates our Monte Carlo random walk ap-proach. To illustrate the above, we show in Fig. 5 simula-tion results for random walkers with both unconstrainedand constrained mobility, generated with step lengths N i n f e c t ed l = 5 〈 r 〉, exponential growth
10 100
Number of time steps (t) N i n f e c t ed l = 〈 r 〉 , power-law growth (A)(B) A = 2498 ± 1B = 1303 ± 57 µ = 0.0646 ± 0.0003 m = 1.87 ± 0.01 A = 2526 ± 4B = 292 ± 43 µ = 0.0377 ± 0.0006m = 2.30 ± 0.05 FIG. 5. Monte Carlo results for (A) exponential and (B) con-strained growth (both with saturation) fitted to a generalizedlogistic function N ( t ) = A (1 + Be − µt ) / (1 − m ) . The data arethe same as shown in Fig. 2 l = 5 h r i and h r i respectively. While the generalized lo-gistic function provides a reasonable fit for the uncon-strained curve (exponential growth), a large discrepancyis observed in the other (power-law) case, with signifi-cantly different values for the fit parameters. It is furtherevident that for t α type power-law growth, the infectionrate dN/dt should be proportional to αt α − . This is alsosupported by data. As an example, we show data cor-responding to three-day averaged values for the reporteddaily deaths from India and their corresponding power-law fit in the top panel of Fig. 6. As expected, we obtain agrowth exponent of β = α − α = 2 . dN/dt following a t α − power-law increase. This observed consistency furtheraffirms the validity of our Monte Carlo method. Thus,our general observations suggest that the growth curvesfrom effectively contained scenarios always ought to befitted accordingly, by including power-law behavior. Thissupports the contention that constrained growth curvesfrom global COVID-19 data necessarily require epidemio-logical analyses that incorporate additional mechanisms,similar to those described in Refs. [6, 8]. Summary and conclusions
In summary, we used asimple two-dimensional random walk Monte Carlo modelto study the spread of Covid-19-like infection within acontained population. Apart from proximity based con-tact, our model has no underlying assumptions about the
10 20 30
Days 〈 d N / d t 〉 da ys Averaged daily deaths (India)
10 20 30 40
Time steps 〈 d N / d t 〉 t i m e s t ep s Simulated infections (6k case) β = 1.82 ± 0.08 β = 1.26 ± 0.09 FIG. 6. Top panel: Growth in three-day-averaged dailydeaths reported from India and its corresponding power lawfit h dN/dt i = At β . Bottom panel: Similar data obtainedfrom our random walk Monte Carlo simulations, for a citywith population of 6k and density 10k/unit area. In bothplots the β are quoted with ± σ uncertainty. nature of infection spread or its reproduction number, re-covery time etc. Three growth regimes, corresponding todifferent levels of containment are shown to emerge nat-urally from our simulations. Under stringent conditions,so that only nearest-neighbor connections are allowed,our simulation results show a power-law growth in time,with growth exponents α = 2 . − .
3, similar to initialCOVID-19 data from China [6]. Furthermore, the growthappears to be consistent, with no apparent dependenceon population size or density. Based on available data,this analysis suggests that the containment and mitiga-tion strategies employed/followed in Chinese provincesafter the initial outbreak were near optimal and resultedin growth exponents that were close to the smallest limit-ing value. On comparison with data from other countries,we observe that reasonable agreement can be attained byintroducing small-world-type connections in the simula-tion model. We anticipate that such a Monte Carlo ap-proach (and its more generalized versions) will be usefulfor the evaluation of future strategies in the midst of thepresent pandemic. As concluding remarks, we briefly mention the generalsimilarity between (i) peripheral growth, that requiresan extra spatial diffusion term [8] in the SIR model, (ii)our simulation results for short step-lengths taken by therandom walkers, and (iii) the diffusion of particles to dis-tinct sites on a two-dimensional lattice [33] at short time-scales. All these cases show a quadratic growth in time.We further observe that a simple logarithmic correctionto our quadratic results (so that t → t ln t ) yields apower-law exponent of about 2.5, in rough agreementwith the intermediate values for contained growth, bothdescribed here and observed in Refs. [6, 7]. Further in-vestigations into this potential connection present an in-teresting research problem for both epidemiologists andphysicists alike. Acknowledgments
We are grateful to Prof. N. Barikfor fruitful discussions and to Prof. S. M. Bhattacharjeefor directing us to Ref. [33] and suggesting the possibilityof logarithmic corrections to t . ST acknowledges fundingsupport from the National Research Foundation, SouthAfrica, under Grant No. 85100. [1] J. T. Wu, K. Leung, and G. M. Leung,Lancet (London, England) , 689 (2020).[2] M. Chinazzi, J. T. Davis, M. Ajelli, C. Gioan-nini, M. Litvinova, S. Merler, A. Pastore y Pi-ontti, K. Mu, L. Rossi, K. Sun, C. Viboud,X. Xiong, H. Yu, M. E. Halloran, I. M. Longini,and A. Vespignani, Science , 395 (2020),https://science.sciencemag.org/content/368/6489/395.full.pdf.[3] B. Tang, F. Xia, S. Tang, N. L. Bragazzi,Q. Li, X. Sun, J. Liang, Y. Xiao, and J. Wu,International Journal of Infectious Diseases , 288 (2020).[4] S. Zhao, Q. Lin, J. Ran, S. S. Musa, G. Yang, W. Wang,Y. Lou, D. Gao, L. Yang, D. He, and M. H. Wang,International Journal of Infectious Diseases , 214 (2020).[5] T. Zhou, Q. Liu, Z. Yang, J. Liao,K. Yang, W. Bai, X. Lu, and W. Zhang,Journal of Evidence-Based Medicine , 3 (2020),https://onlinelibrary.wiley.com/doi/pdf/10.1111/jebm.12376.[6] B. F. Maier and D. Brockmann, Science [11] N. Gan, “Beijing extends residential lockdowns and tight-ens outbound travel as coronavirus infections spread.https://edition.cnn.com/2020/06/16/asia/coronavirus-bejing-spread-intl-hnk/index.html,” (2020).[12] “World Health Organization. https://covid19.who.int/,”(2020).[13] S. Meyer and L. Held, Ann. Appl. Stat. , 1461 (2020).[22] D. J. Watts and S. H. Strogatz, Nature , 440 (1998).[23] A. L. Ziff and R. M. Ziff, “Fractal kinetics ofCOVID-19 pandemic. medRxiv: 2020.02.16.20023820. https://doi.org/10.1101/2020.02.16.20023820,” (2020).[24] O. Hallatschek and D. S. Fisher,Proceedings of the National Academy of Sciences , 290 (1959),https://academic.oup.com/jxb/article-pdf/10/2/290/1411755/10-2-290.pdf.[33] H. Larralde, P. Trunfio, S. Havlin, H. E. Stanley, andG. H. Weiss, Phys. Rev. A45