DDiffusion as a First Model of Spread of Viral Infection
Paulo H. Acioli ∗ Department of Physics and Astronomy,Northeastern Illinois University, Chicago, IL 60625 (Dated: May 27, 2020)
Abstract
The appearance of the coronavirus (COVID-19) in late 2019 has dominated the news in thelast few months as it developed into a pandemic. In many mathematics and physics classrooms,instructors are using the time series of the number of cases to show exponential growth of theinfection. In this manuscript we propose a simple diffusion process as the mode of spreadinginfections. This model is less sophisticated than other models in the literature, but it can capturethe exponential growth and it can explain it in terms of mobility (diffusion constant), populationdensity, and probability of transmission. Students can change the parameters and determine thegrowth rate and predict the total number of cases as a function of time. Students are also giventhe opportunity to add other factors that are not considered in the simple diffusion model. a r X i v : . [ phy s i c s . e d - ph ] M a y . INTRODUCTION The end of 2019 and beginning of 2020 have been dominated by the spread of the coron-avirus (COVID-19). The disease started in the Hubei province in China in December 2019and by early January 2020 started to spread. It grew very rapidly, triggering responses fromthe Chinese and other governments. On March 11, 2020 the World Health organizationdeclared that COVID-19 was then characterized as a pandemic. The situation triggereddifferent responses from different governments. In The United States, many colleges tookthe initiative to start their own social distance programs, including sending all studentshome, extending spring breaks, and ultimately moving all classes to distance learning.
Many cities and states followed up with shelter at home mandates. The situation in manycountries became alarming due to the exponential growth of new cases and deaths.
One of the consequences of this pandemic is that instructors at colleges and universitiesstarted to monitor and model the data, using this as a teachable moment for students andcolleagues. The first step in modeling the data is to plot the number of cases as a functionof time and show that it exhibits an exponential growth. For beginners it is important tointroduce them to the logarithmic scale, where the plot becomes a straight line and studentscan extract the exponent and realize that the fit to the US data on March 20, 2020 showsthat the number of cases doubled every 2.4 days. Fits to the initial data for Chicago andNew York City show that New York has a higher growth rate than Chicago FIG. 1, and thetotal number of cases in Chicago is substantially smaller. One can speculate that this mightbe due to difference in population densities in these cities. After the first 20 days and aftermitigation and social distance measures were imposed, it is clear that the rate of infectionslowed down.In this work we propose to look at some of the factors affecting the spread of viruses usinga simple diffusion model in which each individual in a population is treated as a Brownian particle with diffusion constant D . Also added to this model is the incubation period of thevirus and a probability of transmission of the virus if individuals are closer than a certaindistance. This model is to be used as a project in a computational physics course and verifyif it adequately predicts the exponential growth of the number of cases, as well as if thepopulation density, mobility, and probability of transmission play roles on the percentage ofthe population that will be infected as a function of time. Students, will be asked to modify2 IG. 1. COVID-19 cases in Chicago (CHI,blue) and New York City (NYC,red). The dashedlines are the best fit to the date after the initial spike in cases. The x -axis is the number of daysfrom the date that the first case in each city was detected. a code, analyze the output for different sets of parameters, and write a critical analysis ofits predictive effectiveness. II. COMPUTATIONAL PROJECT DETAILS
In this section we propose a project to be implemented in computational physics courses tostudy the spread of infectious diseases as a simple diffusion of individuals. The benefits of thisproject is that its implementation is simple, but it can lead to a qualitative understandingon how diseases are spread and it can also allow the student to understand factors that canaffect it. The main ingredients of the model are: individuals are considered particles thatobey a Brownian diffusion process; each individual will have three possible states, healthy,sick (contagious), and cured; a healthy individual has a probability of getting infected if itsdistance to a sick individual is smaller than a certain threshold; the incubation and sickness3eriods are the same; once an individual gets cured, it cannot be infected or contagious again.VPython (or GlowScript) is the language of choice as it allows for a real time visualizationof the infection spread and it will allow for fast simulations of small populations on a laptop.In the assignment, students will implement the code, analyze the data generated, critiquethe initial assumptions, and propose improvements for more realistic simulations. Below wedescribe the algorithm starting with the standard diffusion equation.The diffusion equation in 2D ∂f ( x, y, t ) ∂t = D (cid:20) ∂ f ( x, y, t ) ∂x + ∂ f ( x, y, t ) ∂y (cid:21) (1)is used to study many phenomena from diffusion inside the nucleus to population dy-namics to solving the Shr¨odinger equation. The function f will have different meaningsdepending on the application, from temperature to density. In this work we will use it forthe diffusion of individuals, treated as particles, over closed boundaries subject to contami-nation of a viral infection. The normalized solution to Eq. (1) for a single particle is givenby f ( x, y, t ) = f √ πDt e − x y Dt . (2)Therefore one can simulate the diffusion of a particle from its previous position by generatinga Gaussian distribution of zero mean and variance √ Dt . For a system on N non-interactingparticles with the same diffusion constant we use Eq. (2) for each particle at each simulationtime step.The next ingredients in the simulation will be the population density ρ , the numberof habitants (particles in the simulation cell) N pop , the diffusion constant D , the numberof simulation steps N step , the time step dt , the incubation period ( t inc ), the transmissionradius, and the probability of transmission from an infected to a healthy individual ( prob ).All of these variables are set at the beginning of the simulation. We preset that the totalsimulation time consists of 90 days and that each time step is 0.01 days, therefore eachsimulation takes 9000 steps. The algorithm is described below:1. Input N pop , N step , D , dt , t inc , prob , r transm
2. Calculate the size of the square cell as L = (cid:112) ( N pop/ρ ))3. Initialize the initial population 4. Choose a fraction of the initial population to be infected, and set timer for the sickness( t sick )5. Loop over N step
6. Move all individuals according to the Gaussian distribution (Eq. (2)).7. Compute the distance between each healthy and infected individualsIf the distance is less than r transm , the healthy individual becomes sick with prob-ability prob
8. Subtract the sickness timer by dt
9. If t sick < III. AN EXAMPLE OF A PROJECT
In this section we provide an example of a project that students could pursue and alsovalidate the method by comparing the results of the simulation with the data presented inFIG. 1. In this project, students will study the differences of the infection proliferation inNew Your City (NYC) and in Chicago. New York City was chosen because it is experiencinga very rapid growth on the number of cases and preventive measures such as shelter at homewere taken by the city and state governments at very early stages of the infection. Chicagowas chosen as a local connection to Northeastern Illinois University (NEIU) students, and isalso experiencing an exponential growth in the number of cases and has also been affected5y shelter at home mandates. We must be mindful that the current model will not be ableto study the effects of the measures that the government are taking, but it might be able tojustify their need.The population density of New York city is 10,194 people/km , while Chicago’s is 4,665people/km . It is natural to offer the students the hypothesis that if all other variables arethe same, the spread of infections in New York City will be faster than in Chicago and thatthe number of cases will be much larger for the same period of time after the first case. Withthis hypothesis alone, students should be able to generate a enough data for the project.They can also discuss if the mobility (diffusion constant) should be the same in both casesand study the effect of mobility in the spread of the disease. In the proposed model theperiod of incubation is the same as the period of sickness, students can be offered the optionto modify this assumption. The model also assumes that once cured, an individual willattain immunity, and will not be able to spread the disease. They can discuss modificationsto the model to incorporate relapse. It is clear that with this very simple model and thisvery limited two-city project, they can perform a very thorough study that can give us somequalitative understanding of infectious disease spread.In order to validate our model we will try to recreate the same initial exponential growthas experience by Chicago and New York City. The exponential fit to the data in FIG. 1 yieldsexponents of 0.352 and 0.455 for CHI and NYC, respectively. Although we envision that wecan see such growth in a 100 individual population, this growth is unattainable for the sametime frame as the real data, as 100% of the population would be infected in days 13 and 11 ifwe start with one individual infected in day 0. Deviations from the exponential behavior willbe seen even earlier as there is not enough susceptible population available to be infected.In order to get similar exponential growth we start with simulations of populations of 100individuals and 1 sick individual chosen at random. The initial diffusion constant was chosenas D = 100m /day with a time step of dt = 0 . √
2. In this case the number of simulation steps is 9000 for each individualsimulation is 27s in a laptop with Intel(R) Core(TM) i5-8250U CPU @ 1.60GHz. Therefore50 simulations will take about 22.5 minutes of computation. We assume that the radius ofcontamination is 2m, and that the probability of contamination is 20% per time step. Theincubation period is taken as 14 days. Below are some of the results of these simulations.In FIG. 2a we present the results of the simulation for New York City with the set of6arameters in the previous paragraph. As one can see, under this assumption 50% of thepopulation is sick on day 11 and 99% of the population will be infected after 28 days. Onecan see that after day 18 the number of sick people starts to decline. However, if the deathrate is similar to what has been observed on the COVID-19 pandemic, about 5.5% of thepopulation of NYC would perish, and the numbers would be even worse, since no major cityin the world would be able to have hospital beds for 83.8% of its population at the peakof the infection. In the figure we also observe that the rate of cure follows the number ofinfected individuals with a lag time of 14 days, which is the incubation/sickness period. Anexponential fit, shown in FIG. 2b, yields an exponent of 0.41 which is a bit lower than whatwas observed in the real data. Students will need to search for the set of parameters that isable to fully match the initial growth rate data that they choose to model.In FIG. 3a we report the results for the city of Chicago. Most of the parameters arethe same, with the exception of the density that will be changed to ρ = 0 . .For this set of parameters 50 % of the population will be infected after 26 days reaching amaximum of 87.8% after day 86. The peak of the number of sick people is 42.3 % and ithappens about 32 days after the first 1% of the population is infected. In FIG. 3b we presenta logarithmic scale plot of the total number of infections together with and exponential fitfor the initial stages of the infection. Similarly to what happened in the simulation for NYCit is clear that the number of infections deviates from the exponential growth after 13 days.The exponent in the fit is 0.2 which is also smaller than the value for the data for Chicago.A closer fit is likely to take place by increasing the probability of infection.These two examples do support our initial hypothesis that under the same conditions, onewould expect a slower growth rate in a less dense population. This leads to the conclusionthat keeping all parameters the same, the percentage of the population infected is correlatedto the population density. In order to determine how mobility affects the rate of infection, onewould change the diffusion constant and repeat the simulations. The expectation is that thelower the mobility, the lower the total infection rate. In addition, even with the limitationsof a small population simulation we were able to find growth rates of infections that aresimilar to what was observed for the cities of Chicago and New York. Thus validating thisapproach to simulate the spread of infectious diseases.In order to show the limitations of small population simulations we run a simulationfor the city of Chicago with a probability of infection of 30% per time step but increased7 IG. 2. a) Average of 50 90-day simulations of the spread of a virus in a population of 100individuals in a square cell with the population density of NYC ( ρ = 0 .
012 people/m ). D =100m /day, prob = 0 . dt = 0 . the number of individuals to 10000. The results are shown in FIG. 4a. With this set ofparameters we obtain the same exponential growth as observed in FIG. 1 for a period of15 days. Starting with 0.01 % of the population infected at the end of 90 days about 40 %of the population will have been infected. One can see that even a small diffusion constantof 100 m /day and a probability of transmission of 30 %, the virus will spread to a largepercentage of the population. It is clear that a diffusion model does a very good job ofsimulating the spread of the COVID-19 virus and could be used in a classroom setting, andwith the modifications suggested below one can even use it for more realistic predictions when8 IG. 3. Average of 50 90-day simulations of the spread of a virus in a population of 100 individualsin a square cell with the population density of the city of Chicago ( ρ = 0 . ). D =100m /day, prob = 0 . dt = 0 . simulating larger populations. The main limitation will be access to better computationalresources such as a parallelized code and a computer cluster.To finalize, we show snapshots of the simulation cell and its population in FIG. 5. Whitespheres represent the susceptible people, red represent sick, and green are those individualsthat recovered. This window is very useful for a quick analysis of what is happening andcan show how the infection is spread over time. We found it useful to change the radius ofthe individuals so that they are visible as the population increases, and the cell size on thescreen remains the same. Because we are displaying the results for 10000 individuals it is9 IG. 4. One 90-day simulations of the spread of a virus in a population of 10000 individuals ina square cell with the population density of the city of Chicago ( ρ = 0 . ). D =100m /day, prob = 0 . dt = 0 . interesting to see that the infection has a clear origin and the diffusion process spreads thedisease outwards. In these simulations we chose closed boundaries and as a result there willbe a limitation on the growth as there will be less susceptible neighbors to be infected whenthe disease reaches a boundary. To avoid this limitation one can use periodic boundaryconditions.The main limitation of this model is that it assumes the same mobility for all individualsand a constant density in the simulation cell, and no travel between cities. In addition,the average distance traveled by Brownian particles is proportional to the square of the10 IG. 5. Snapshots of one of the simulation 90-day simulations of the spread of a virus in apopulation of 10000 individuals in a square cell with the population density of CHI ( ρ = 0 . ). D = 100m /day, prob = 0 . dt = 0 . simulation time. Thus, there is a limit on how far they will travel and how many peoplean individual can infect. However, depending on the level of the students, instructors canuse the same code and create sub-regions in the simulation cell that would have differentpopulation densities and subsets of the populations can be confined to these regions. Thestudents could allow individuals to go from one region to another with a given probability.This would allow the possibility of an infected individual to move to a region that otherwisewould have no infections. This modification would allow for the spread of the infectionacross borders. One could also include quarantine effects, by creating small cells where asingle infected individual will be confined for a period of time and no other individuals areallowed in. Social distancing can be implemented by a small repulsive potential attributed tosome individuals. These modifications will bring this simple model closer to more realistic11imulations as those of references . To help the implementation of these projects, theoriginal code can be downloaded from the author’s website . IV. CONCLUSION
This project on itself is simple enough that it can be implemented and analyzed in afirst course on computational physics. However, it is rich enough that can lead to a largeamount of data that can be used to qualitatively and quantitatively analyze the spread ofviral infection. In addition, it can help students understand the phenomenon of particlediffusion and Brownian motion. In the example presented, the simulation was limited to thesame 20% chance of infection, the same time step, the size of the population was fixed, andthe length of the disease was the same as the incubation period. Under these conditionsstudents will have a larger number of parameters to change and verify if they can reproducedata readily available in the news. In addition, they will have the opportunity to determinethat the initial growth rate in the number of cases follows an exponential trend and noticethat once it reaches a certain threshold if it will follow the expected logistic behavior.In conclusion, we hope to have convinced the reader that a simple diffusion model canbe used to qualitatively explain the spread of disease and even in some cases quantify it. Itis our judgment that it will allow students to work on a problem that is directly affectingthem and that these simulations will help them offer valuable insight on the issue.
ACKNOWLEDGMENTS
I would like to thank Dr. Greg Anderson for sparking the discussion of using the COVID-19 pandemic as a teachable moment and Dr. Orin Harris for inspiring the idea of comparingthe data for different populations. ∗ [email protected] https://coronavirus.jhu.edu/map.html Feynman, R. (1964). ”The Brownian Movement”. The Feynman Lectures of Physics, Volume I.pp. 411. https://vpython.org/ F. A. Williams, Amer. J. Phys. , 467 (1958). D. C. Kelly, Amer. J. Phys. , 585 (1968). Karmeshu, and L. S. Kothari, Amer. J. Phys. , 1264 (1972). C. Domb and E. L. Offenbacher, Amer. J. Phys. , 49 (1978). J. DArruda, and E. W. Larsen, Amer. J. Phys. , 392 (1978). A. M. Albano, N. B. Abraham, D. E. Chyba, and M. Martelli, Amer. J. Phys. , 161 (1984). S. Redner and P. L. Krapivsky, Amer. J. Phys. , 1277 (1999). K. Ghosh, K. A. Dill, M. M. Inamdar, E. Seitaridou, and R. Phillips, Amer. J. Phys. , 123(2006). D. R. Spiegel and S. Tuli, Amer. J. Phys. , 747 (2011). T. Pang, Amer. J. Phys. , 980 (2014). G. Rossetti, L. Milli, S. Rinzivillo, A. Sˆırbu, D. Pedreschi, and F. Giannotti Int. J. Data Sci.and Anal. , 61(2018). https://towardsdatascience.com/modelling-the-coronavirus-epidemic-spreading-in-a-city-with-python-babd14d82fa2 The code can be downloaded from http://physics.neiu.edu/research/viraldiffusion/Viral Diffusion.py. nstructors and students using the code should cite the present manuscript in any presentationor publication that uses this code as downloaded or modified.nstructors and students using the code should cite the present manuscript in any presentationor publication that uses this code as downloaded or modified.