Responding to COVID-19 under uncertainty: a simulation based analysis

I became concerned by a discussion I heard related to COVID-19 this past weekend. The debate centered on whether or not it was a good idea to respond to COVID-19 progressively as new information became available, or if it was better to first respond energetically regardless of what we know at this point about the pandemic. I was surprised by the discussion because people arguing about it were building their arguments based on the premise that they knew exactly what was going on. This is not a surprise as health officials (i.e. at least in Mexico) are building the case for their response -in this case waiting a bit more to enforce some serious social distancing measures- based on the premise that they know how much people is being infected and how infectious and deadly this virus is. Yet, we simply do not know any of these parameters, and betting on we knowing could potentially lead to some serious mistakes.

It is possible to use mathematical modeling and numerical simulation to look more deeply into these arguments. For this, let’s consider a pandemic outbreak in a population of 120 million people which behaves somewhat similar to COVID-19 in terms of how quickly it spreads and in its -so far- known mortality rate. Further, let’s also assume it is possible to model this pandemic using a slightly modified version of the classic SIR (i.e. S: Disease, I: Infectious and R: Recovered) model. The two additions I have made to this model are: 1) there are information delays between the real number of infected people and the confirmed number of cases of the disease, 2) the mortality rate of the disease increases as it overburdens the health system. For the latter point I have chosen a very conservative effect of the overburden of diseased patients on hospitals' effectiveness. This is: if the number of infected people overpasses the carrying capacity of the health system (in this case 100 thousand patients), then the mortality rate of this pandemic increases by twenty percent. This is of course not the case, it could be much higher if overcrowding of hospitals becomes very severe. But, let’s shy away from alarmism, this is not the point of the analysis. If you want to have a look at the model, by all means do so, the model is at the bottom of this blog.

Let’s first consider the effect of social distancing. The model reproduces the effect that social distancing has on flattening the curve of infected people, as shown in the graph below:

Three cases of social distancing are being displayed here. None (blue), indicates that people can go on with their lives without making any changes, mild (orange), that they do implement some adjustments (i.e. they reduce their contact rate with other people by 30%) and aggressive (red), that they implement some serious noticeable changes (i.e. they reduce their contact rate with other people by 50%). As you can see, the more aggressive we are with social distancing, the longer it takes for the virus to spread through the population, potentially preventing the overcrowding of hospitals.

Thus social distancing is key to respond to this outbreak. However, when is it a good time to implement this policy and how relevant is the fact that we simply do not know how many people are infected at any point in time? It turns out this is very important. The graph below highlights this for a case in which it takes on average 15 days to track suspected cases, apply the tests and confirm who is infected and who is not.

The dark blue line tracks the real number of people that has been infected with the disease, the lighter blue line tracks the number of confirmed cases. This is the number that officials are using for reporting on the situation. Please note the difference at day 120 of the outbreak. At this point, considering a 15 day average delay time, the difference between the confirmed number of cases and the real cases is 273%.

A lot of factors determine this average delay time, including the level of sophistication of the nations’ health infrastructure, the readiness of emergency response agencies and, I believe very importantly, the number of tests being applied to the population. If you are not testing the population, it will be very difficult to build an accurate picture of what is going on. Consider, for instance, the simulations below:

Each panel shows the difference between the number of confirmed and real cases for three different information delay times. The left panel (1 day delay time) is the ideal case in which you are able to track the progression of the disease with great detail, perhaps this is close to what happened in South Korea and Taiwan. The panel to the right displays a terrible a case, in which it takes an entire month, figuring out what is going on. In this case, at day 120, the difference between the confirmed and real cases is 480%.

Reducing this information delay through testing and preparedness is incredibly important from a policy perspective, as it determines when it is necessary to trigger containment actions (e.g., social distancing). For seeing this, let’s assume that containment actions are triggered once the number of confirmed cases surpasses 100 thousand cases (0.10 million) and let’s zoom in to the first days of the outbreak as shown below:

The color legend this time denotes different information delay times. If the response to the outbreak is running a 1 day information delay (blue), containment actions are triggered at day 86, for 15 days (orange) at day 94 and for 31 days (red) at 98. In the latter case, you will be running 12 days behind the disease.

This is important because 12 days have passed without implementing potentially lifesaving actions. The graph below shows just how important social distancing can be to save lives.

The color legend denotes the three social distancing approaches considered. Note that both mild and none social distancing increase the lethality of the disease, while aggressive social distancing in fact reduces its lethality (right axis). The left axis denotes the actual number of people diseased. The stark contrast of the three lines highlights what is at stake here. However, do not take these numbers at face value, this is an exploratory analysis, not a prediction. No simulation analysis should be used as such. Below you can explore in more detail the assumptions and limitations behind these estimations.

We can do a similar exercise to highlight the importance of information delays in responding to the outbreak as shown below:

The trend is quite clear: the longer it takes you to collect the necessary information to mount an appropriate response, the deadlier the disease becomes. Consider once more that I am modeling quite conservatively the effect of overcrowding on the health system. If I am terrible underestimating this effect, this line becomes stepper.

The preceding discussion provides the necessary background for the actual policy discussion. This is: what is the difference between responding sequentially or in one step to this outbreak? For this I capture the first policy-sequential policy- by assuming that social distancing is implemented in two phases. In the first phase, only half the effort of social distancing is implemented and in the second phase, the full effort of social distancing is implemented. Thus, if the decision is to implement aggressive social distancing. Under this policy, during the first phase, mild social distancing is implemented, and once the conditions of the second phase are achieved (measured by the confirmed number of cases), then aggressive social distancing is implemented. The first phase is triggered when 10,000 cases are confirmed. The second phase is triggered when 100,000 cases are confirmed. The second policy –one step policy- is modeled in a much simpler way: social distancing is implemented when 10,000 cases have been confirmed. The following panel shows how these policies behave across two metrics: a) mortality rate and b) number of days with social distancing in effect (i.e. this is proxy for economic fallout), for two cases of social distancing (i.e., aggressive social distancing and mild social distancing).

The square (blue) dot tracks the performance of the sequential policy, the circle (orange) dot shows the performance of the one-step policy. If aggressive social distancing is to be implemented (top panel), then there is no difference in the number of days required for containment between both policies. Yet, the one-step policy results in a lower mortality rate (saving approximately 8,000 more lives). This effect is accentuated with longer information delays. Thus, nations that are less ready to face the outbreak and that are seeking to implement aggressive social distancing are better served by doing this in one-step. If, instead, mild social distancing is implemented, both policies result in a mortality rate higher than two percent, albeit the one-step policy still achieves a slight lower mortality rate (i.e., approximately 1,000 less casualties). However, in this case, the one-step policy results in a longer time required for containment (i.e. approximately 2 more weeks). Nations seeking to implement mild social distancing to respond to the outbreak then need to balance this trade-off. If you ask me, saving one thousand people justifies without question having to put up with a bit more time of extra mild social distancing. Thus in my view, so far there is no strong evidence that supports favoring the sequential policy over the one-step policy.

Yet, these results are of course dependent on various uncertain parameters. Two seem to be particularly relevant and uncertain right now: 1) the infectivity and 2) the mortality rate of COVID-19. There is some talk that some climates are less favorable to the spread of COVID-19 and the data coming from Asia, Europe and the US it is still inconclusive with respect to how deadly this novel virus is. Then, the question is whether or not conclusion above still holds when considering different scenarios of this outbreak.

We can explore this argument by running the model a few thousand times across variations in infectivity, death rate, average delay time and social distancing. The results of this experiment are shown below:

The color legend in this heatmap denotes the population fatality rate for this outbreak across scenarios that combine different virus infectivity and mortality rate values. The percent value of this parameter indicates how much worse (above 100%) or more favorable (below 100%) the virus turns out to be with respect to the baseline parameters analyzed above. The four panels combine policies and social distancing implementation. The left columns show results for the sequential policy, the right columns show results for the one-step policy. The top row shows results for aggressive social distancing, the bottom row shows results for low social distancing.

If we focus only on how green each panel is, then it is clear that the one-step policy with aggressive social distancing outperforms the other three combinations. This experiment also shows that implementing either policy with mild social distancing is not wise as it exposes the population to higher risks.

The results also highlight the impact of contextual factors in this discussion. Notice that when the infectivity rate is half of what we originally estimated, this is: only 5% of contacts between infected and uninfected result in contagion, the fatality of the outbreak is contained regardless of policy choice. In others words, if this were to happen, major disruptions would be avoided, not due to a talented policy intervention, but to mere luck. Thus, results of this crisis should be weighted considering not only what was the end result, but also considering how those results came to be.

This brief and basic analysis shows that under uncertainty it is more robust to take social distancing and testing seriously and to be aggressive about these two policy elements, even as the crisis is beginning to unfold. In this way, we do not depend so much on luck to avoid irreparable lost.

You can use the dashboard below to interact with the analysis:

The code for this this analysis can be found here:

covid.epidemic <- function(t, state, parameters) { with(as.list(c(state,parameters)), { #Endogenous auxiliary variables if (policy==1) { Social.distancing<-ifelse( PolicySwitch_1>0, ifelse(PolicySwitch_2>0,Social.distancing.policy,min(Social.distancing.policy*1.6,1.0)) ,1.0 ) } else if (policy==2) { Social.distancing<-ifelse(PolicySwitch_1>0,Social.distancing.policy,1.0) } else if (policy==3) { Social.distancing<-Social.distancing.policy } else { Social.distancing<-1.0 } Probability.of.Contact.with.Infected.Person<-population.infected.with.COVID/Total.Population Susceptible.Contacts<-population.susceptible.to.COVID*Contact.Frequency*Social.distancing Contacts.between.Infected.and.Uninfected.People<-Susceptible.Contacts*Probability.of.Contact.with.Infected.Person #[people] #Flow variables Infection.Rate<-Contacts.between.Infected.and.Uninfected.People*Infectivity #[people/day] Recovery.Rate<-(population.infected.with.COVID*(1-mortality.rate))/Average.Duration.Of.Infectivity population.infected.with.COVID.that.requires.hospitalization<-population.infected.with.COVID*hospitalization.rate health.system.load<-population.infected.with.COVID.that.requires.hospitalization/Health.system.carrying.capacity effect.of.overburden.on.mortality.rate<-ifelse(health.system.load>1,1.2,1.0) Death.Rate<-(population.infected.with.COVID*mortality.rate*effect.of.overburden.on.mortality.rate)/Average.Duration.Of.Infectivity #Information delay dperceived.population.infected.with.COVID<-(population.infected.with.COVID-perceived.population.infected.with.COVID)/average.delay.time #State variables dPolicySwitch_1<-ifelse(perceived.population.infected.with.COVID > Social.distancing.trigger.s1,1,0) dPolicySwitch_2<-ifelse(perceived.population.infected.with.COVID > Social.distancing.trigger.s2,1,0) ddeaths<-Death.Rate dpopulation.susceptible.to.COVID<-(-1)*Infection.Rate # [people] dpopulation.infected.with.COVID<-Infection.Rate-Recovery.Rate # [people] dpopulation.recovered.from.COVID<-Recovery.Rate # [people] list(c(dpopulation.susceptible.to.COVID, dpopulation.infected.with.COVID, dperceived.population.infected.with.COVID, dpopulation.recovered.from.COVID, ddeaths, dPolicySwitch_1, dPolicySwitch_2)) }) } InitialConditions <- c( population.susceptible.to.COVID = 120 , population.infected.with.COVID = 1/1e6, perceived.population.infected.with.COVID = 0, population.recovered.from.COVID = 0, deaths=0, PolicySwitch_1=0, PolicySwitch_2=0 ) times <- seq(0 , #inicial time #days 500 , #end time #days 1 ) #time step #days intg.method<-c("rk4") #==================================================================== #Experimental design Part 1 #==================================================================== Exp.design<-expand.grid( seq(0.1,6.0,by=0.2), c(0.5,0.7,1.0) ) colnames(Exp.design)<-c("average.delay.time","Social.distancing.policy") Exp.design$Future.ID<-c(1:nrow(Exp.design)) #add policy tag Exp.design$Policy<-0 #Duplicate for rest of policies #policy 1 Exp.design_p2<-Exp.design Exp.design_p2$Policy<-1 #policy 2 Exp.design_p3<-Exp.design Exp.design_p3$Policy<-2 #Rbind all Exp.design<-rbind( Exp.design, Exp.design_p2, Exp.design_p3 ) #add run ID identifiers Exp.design$Run.ID<-c(1:nrow(Exp.design)) #execute runs runs<-apply(Exp.design,1, function(x) { parameters<-c( Infectivity = 0.1, # [1] dimmensionless mortality.rate = 0.02, # [1] dimmensionless Contact.Frequency = 2, # people/day Total.Population = 120, # million people, average.delay.time = 10*as.numeric(x["average.delay.time"]), #days Social.distancing.trigger.s1 = 0.01, #million people Social.distancing.trigger.s2 = 0.1, #million people Social.distancing.policy = as.numeric(x['Social.distancing.policy']), Health.system.carrying.capacity = 0.1, #million people Average.Duration.Of.Infectivity = 15, # days hospitalization.rate = 0.2, #[1] dimmensionless policy=as.numeric(x['Policy'])); #days out <- ode(y = InitialConditions, times = times, func = covid.epidemic, parms = parameters, method =intg.method ); out<-data.frame(out); out$Run.ID<-x['Run.ID']; out } ) # runs<-do.call("rbind",runs) #merge runs with experimental design dim(runs) runs<-merge(runs,Exp.design, by = "Run.ID",all.x=TRUE) dim(runs) # dir.data <- "~\\Edmundo-ITESM\\2.Cursos Impartidos\\Modelacion de Sistemas\\2020\\"

write.csv(runs,paste0(dir.data,"covid.csv",sep=""),row.names=FALSE) #now aggregate results #mortality of pandemic sruns<-subset(runs,time==500)[,c('Run.ID','deaths')] #number of days until normality runs$days.out<-ifelse(runs$population.infected.with.COVID>0.01,1,0) sdays<-aggregate(runs$days.out,list( Run.ID=runs$Run.ID ), sum ) sdays$days.out<-sdays$x sdays$x<-NULL #merge both summary statistics dim(sruns) sruns<-merge(sruns,sdays, by = "Run.ID") dim(sruns) #merge with experimental design dim(sruns) sruns<-merge(sruns,Exp.design, by = "Run.ID",all.x=TRUE) dim(sruns) # dir.data <- "~\\Edmundo-ITESM\\2.Cursos Impartidos\\Modelacion de Sistemas\\2020\\" write.csv(sruns,paste0(dir.data,"scovid.csv",sep=""),row.names=FALSE) #==================================================================== #Experimental design Part 2 #==================================================================== Exp.design<-expand.grid( c(0.1,1.5,3.0), c(0.5,0.7,1.0), seq(0.5,1.5,by=0.10), seq(0.5,1.5,by=0.10) ) colnames(Exp.design)<-c( "average.delay.time", "Social.distancing.policy", "infectivity", "mortality.rate" ) Exp.design$Future.ID<-c(1:nrow(Exp.design)) #add policy tag Exp.design$Policy<-1 #Duplicate for rest of policies #policy 1 Exp.design_p2<-Exp.design Exp.design_p2$Policy<-2 #Rbind all Exp.design<-rbind( Exp.design, Exp.design_p2 ) #add run ID identifiers Exp.design$Run.ID<-c(1:nrow(Exp.design)) #runs runs<-apply(Exp.design,1, function(x) { parameters<-c( Infectivity = 0.1*as.numeric(x['infectivity']), # [1] dimmensionless mortality.rate = 0.02*as.numeric(x['mortality.rate']), # [1] dimmensionless Contact.Frequency = 2, # people/day Total.Population = 120, # million people, average.delay.time = 10*as.numeric(x["average.delay.time"]), #days Social.distancing.trigger.s1 = 0.01, #million people Social.distancing.trigger.s2 = 0.1, #million people Social.distancing.policy = as.numeric(x['Social.distancing.policy']), Health.system.carrying.capacity = 0.1, #million people Average.Duration.Of.Infectivity = 15, # days hospitalization.rate = 0.2, #[1] dimmensionless policy=as.numeric(x['Policy'])); #days out <- ode(y = InitialConditions, times = times, func = covid.epidemic, parms = parameters, method =intg.method ); out<-data.frame(out); out$Run.ID<-x['Run.ID']; out } ) # runs<-do.call("rbind",runs) #merge runs with experimental design dim(runs) runs<-merge(runs,Exp.design, by = "Run.ID",all.x=TRUE) dim(runs) # #mortality of pandemic sruns<-subset(runs,time==500)[,c('Run.ID','deaths')] #number of days until normality runs$days.out<-ifelse(runs$population.infected.with.COVID>0.01,1,0) sdays<-aggregate(runs$days.out,list( Run.ID=runs$Run.ID ), sum ) sdays$days.out<-sdays$x sdays$x<-NULL #merge both summary statistics dim(sruns) sruns<-merge(sruns,sdays, by = "Run.ID") dim(sruns) #merge with experimental design dim(sruns) sruns<-merge(sruns,Exp.design, by = "Run.ID",all.x=TRUE) dim(sruns) # dir.data <- "~\\Edmundo-ITESM\\2.Cursos Impartidos\\Modelacion de Sistemas\\2020\\" write.csv(sruns,paste0(dir.data,"scovid_2.csv",sep=""),row.names=FALSE)