Virus Spread Simulation

So the virus is spreading. At first our responses appeared to be reactive rather than proactive. We seemed to believe that somehow we would avoid disasters happening elsewhere in the world despite doing exactly the same thing. No-one seemed willing to consider what will be happening in a week or two, let alone the long-term end game. Now finally it seems people are at least starting to take this seriously (as of 23/3/2020). But should we be doing more?

People suggested I translate the numbers from my simulations below into figures for Western Australia. I'm loath to do so, as the model has so many assumptions. But they are probably as good estimates as any, given our large uncertainties. So here they are:

  • Act slowly and halfheartedly*: ~130,000 people die in WA, the virus burns out in about six months and we can relax restrictions
  • Act fast and halfheartedly : ~40,000 people die, the virus continues to trickle through the population for more than a year before we can relax restrictions
  • Act slowly but very strongly**, achieve and maintain suppression: ~60,000 people die, we are clear in about six months, but we then have to maintain low numbers until there is a vaccine
  • Act very fast and very strongly, achieve and maintain suppression: ~600 people die, we are clear in a few months, but have to maintain low numbers until there is a vaccine
*'halfheartedly' means strongly enough to 'flatten the curve' but not strongly enough to achieve suppression**very strongly means strongly enough to achieve suppression

Which one would you chose??

The point of these numbers is not to panic people. The biggest number is possible if most of WA get the virus, and we have death rates like those seen in places where health systems have been overwhelmed, as will happen if we don't slow down the curve. In reality we will not let that happen. But the point is that early action and strong action will save many many lives AND (if we get organised) allow us to get back to a more normal life relatively quickly. The nature of exponential spread and long time lags means that even small delays will be very very costly.

For some very good analysis of the COVID19 pandemic and the best ways to deal with it, I highly recommend this article and this article.

For more details on the simulation, read on...

A virus spread simulation I created, with help from my wonderful son Kiet. The image on the left represents a community of 4000 people, moving around and interacting. Each dot is a person. Black dots are people that are uninfected and susceptible. At the beginning everyone is uninfected and susceptible. This community knows that there is a virus epidemic 'somewhere else' in the world. They have implemented checks/testing/quarantine on incoming people, and they are finding some (these are not shown in the simulation). However, eventually someone who is infected but non-symptomatic (latent) slips through the net and starts circulating within the community. At this point they are shown in the cyan colour. They then become infectious (their colour turns to red) and they start infecting other people, who are then latent (cyan) to start with. Since the community thinks they are catching all infected incomers (and/or they don't have enough tests) they are not testing people in the general community, even when they get sick. It's probably just a normal cold after all. The infection spreads unchecked for a while. It's only when a few people get really really sick that they finally start testing and detecting infections. At this point, the health care system is holding fast and there have been no deaths. By now patient zero and a few others have recovered (they are now green).

There are a few important things to notice at the start. These are maybe clearest in the graph on the right, which just shows how the numbers of dots on the left change over time. The number of detected cases (grey) is much lower than the number of actual infectious people (red), and the number of actual infectious people is much lower than the number of latent infected people (cyan). The number of infected people is also rising at an ever increasing rate, which means things very quickly spiral out of control. This combination of lack of widespread testing, time lag due to asymptomatic infections and ever increasing rates of growth means that even when just a few detentions have occurred, big big problems are already very hard to avoid.

Now let's keep the simulation going...

Finally the community realises that there is a real problem and community transmission is occurring. It still takes a few days, but eventually they act and start social distancing and restricting movement and they also finally start testing people who are sick in the general community and quarantining the ones they detect. Notice how the movement of the people slows down in the animation on the left, and the 'response initiated' line appears on the right. The real number of total infections is hugely higher than the number of detected cases at this point. The health system is already under real strain and nearing tipping point. The social distancing and movement restrictions do have an immediate effect (let's assume that people are now getting serious) on the number of new infections and the number of latent cases starts to fall. But this is because they are now becoming infectious and so the number of infectious people continues to rise. Fast. The health system is soon overwhelmed. The mortality rate increases as doctors are forced to choose who to save with the limited resources they have, and the number of deaths rises much faster than at the start.

Eventually (after quite a while) the slowing of new infections translates into a reduction in sick people. The virus is so well established now that it continues to spread, but more slowly. The health system regains some kind of control, but still needs to triage to some extent for a considerable time. When the virus has spread through most of the population it eventually slows further and (in this community) dies out within a year. Some lucky people never got the disease at all (there are still a few black dots at the end, and the total cases doesn't equal 4000). Many people die. The majority recover and live happily ever after (although their share portfolio and their loved ones may not).

Compare this scenario with the scenarios below.

This one is the same as the first except there is more community testing at the beginning and the response to start social distancing and restricting movement is initiated much faster. What differences do you notice? Number of deaths? Case death rate? Note that in this one the simulation ends after a year, but the virus is still spreading.

This next one is the same as the first, with the same start of social distancing and restricting movement, but the social distancing and movement restriction response is much stronger. What differences do you notice? The case death rate is high, but the infection is eradicated with only a portion of the population affected, and so total deaths are lower than the first one. But assuming that the virus is still spreading outside the community, the community here is now still very vulnerable to new infections coming in, in contrast to the other cases where most people get infected and recover by the end (and are thus hopefully at least somewhat immune).

And finally, this one is what happens if there is good community-wide testing, the response is activated fast (after the second traceable positive test) and the response is strict.

This is just a simulation, but I believe it captures the most important aspects of a new virus invading a community. I tried to set parameter values to represent COVID-19 where I could. This is a small community, which makes it easier to simulate, but you could think of it as a part of a larger city, and multiply the numbers up accordingly. When I did try representing bigger and smaller communities, the results are qualitatively much the same. Widespread community testing (not just new arrivals), and very early action saves many lives. The disease can be suppressed if you go hard enough, but then it will bounce back if you relax the restrictions, unless you totally isolate forever more, or find a vaccine. I would say the best strategies are

1) widespread community testing as much as possible and as early as possible

2) reduce infection rate through strict measures as early as possible (less universal measures only work if there is widespread community testing so you can trace every infected individual, otherwise you probably have to shut everything down)

3a) when the infection rate stabilises at a low level (and not before, and this again depends on widespread community testing to even know) and the number of sick people stabilises at a level that the health system can handle, then consider very slowly and carefully relaxing restrictions, while continuing to test widely anyone with symptoms in order to monitor infection rate and being fully prepared to ramp restrictions back up as soon as there is any hint of things getting out of control


3b) drive for complete eradication by keeping the restrictions on until there is good evidence for zero infection, and then relax restrictions while maintaining high vigilance and levels of testing to detect any new infections as fast as possible (and immediately reapply restrictions if detections occur)

The model is implemented in R and the code is provided below under a GNU GPLv3 license. Please use it and play with the parameters and modify as you wish and see what you discover. If you don't like my assumptions then you are free to change them. To run the code you just need to

1) download the free R software - the basic GUI works best for these kinds of simulations

2) open a new R script from the 'File' menu

3) copy the R code below into the blank script

3) run the script from the 'Edit' menu

### IBM epidemic spread Michael Renton 13/3/2020### GNU GPLv3 General Public License
ws=180n=4000ndaysmax=365latenttime=10recoverytime=16movement=1distance=3pinf=0.16pdeath = 0.00004responsethreshold = 1healthcapacity = 50ptest=0.02dt=0.1

pop = data.frame( x=runif(n)*ws, y=runif(n)*ws, status=0, time=0, detected=FALSE)
#status 0:susceptible, 1:infected latent, 2:infected clinical and infectious, 3: recovered and immune, 4: diedpop$status[1]=1head(pop)absdiff = function(a,b) abs(a-b)

for (t in seq(0,ndaysmax,by=dt)){ ## management if (sum(pop$detected)>responsethreshold & respond==99999999) { movement=0.1 pinf=0.02 ptest=0.2 respond=t } ## move dx = rnorm(n)*movement*(pop$status<4) dy = rnorm(n)*movement*(pop$status<4) dx[pop$detected & pop$status!=3] = 0 dy[pop$detected & pop$status!=3] = 0 pop$x = (pop$x+dx)%%ws pop$y = (pop$y+dy)%%ws ## infect? infected = subset(pop,status==2) whichsus = which(pop$status==0) susceptible = pop[whichsus,] dx=outer(infected$x,susceptible$x,absdiff) dy=outer(infected$y,susceptible$y,absdiff) dd2=dx^2+dy^2 ninf = apply(dd2<distance,2,sum) probinf = 1-(1-pinf)^ninf ifinf = runif(length(probinf ))<probinf pop$status[whichsus[ifinf]]=1 ## change status pop$time[pop$status %in% 1:2]=pop$time[pop$status %in% 1:2]+1 clinical = (pop$status==1 & pop$time>(latenttime/dt)) pop$status[clinical]=2 recover = (pop$status==2 & pop$time> (recoverytime + latenttime)/dt) pop$status[recover]=3 ### mortality nsick=sum(pop$status==2) pdeathnow = pdeath + plogis((nsick-healthcapacity)/20)*15*pdeath pop$status[pop$status==2 & pdeathnow>runif(n)] = 4 ### detection pop$detected[pop$status==2 & runif(n)<ptest]=TRUE pop$detected[pop$status==4 & runif(n)<0.2]=TRUE ## plot if (round(t)==t) png(filename = paste("virusplot",t+1000,".png",sep=""),width = 480*2, height = 480) par(mfrow=c(1,2)) deathrate = round(sum(pop$status==4)/(sum(pop$status>1))*100,2) labb=paste('deaths:',sum(pop$status==4),' case death rate', deathrate ,"%") with(pop,plot(x,y,pch=16,col=c(1,5,2,3,4)[status+1],xlim=c(0,ws),ylim=c(0,ws),xlab=labb)) text(pop$x[1:4],pop$y[1:4],c('patient0','a','b','c'),col=c(1,5,2,3,4)[pop$status+1]) thisrec=c(n-sum(pop$status==0),sum(pop$status==1),sum(pop$status==2),sum(pop$status==3),sum(pop$status==4),sum(pop$detected)) rec=rbind(rec,thisrec) matplot(rec,type='l',col=c(6,5,2,3,4,8),xlab='days',xaxt='n',lwd=2) axis(1,labels=1:ndaysmax,at=(1:ndaysmax)/dt) legend('topleft',c('latent','infectious','recovered','died','total infections','detected infections'),lwd=2,col=c(5,2,3,4,6,8),lty=c(2,3,4,5,1,1)) abline(v=respond/dt,lty=3) text(respond/dt,responsethreshold ,'response initiated') if (round(t)==t) ## finished?? if (sum(pop$status %in% 1:2)==0) break}