1
$\begingroup$

I'm trying to model COVID-19 in New York, and in particular to model the death rate in light of the pre- and post-NY-On-Pause regimes. So I essentially have two SIR models running in parallel: one that reflects the pre (let's call it A) and one that reflects post (call it B). Obviously, A has a higher Beta and B a lower Beta, reflecting little and strict social distancing, respectively.

So I start with the population of New York in A's S, with a few in A's I to get the epidemic started. I then run A and B (which is empty and doing nothing) in parallel until the 23rd day of modeling, at which point I transfer 80% of A's S compartment to B's S compartment, and similarly for the I and R compartments. If I transferred 100% (or 0%), I can see that the combined models would still be modeling what I expect, but I'm a little worried about having two models with substantial populations running in parallel but not interacting after the transfer.

On a practical level, I don't yet know how to use more advanced features of the package I'm using (SimInf in R) to maintain an ongoing interaction, so I did it as best as I could with basic features. On a statistical/modeling level, I've justified this design by saying "Well, people in B are practicing strict social distancing, so wouldn't be likely to interact with people in A who are not." But I worry that there should be a small interaction, which might grow to be significant when A's I grows very rapidly. That is, the number of infectious in A will continue to rise rapidly and will probably exceed B's infectious peak, even with only 20% of the population remaining in A, and that many infectious -- even with low A-B interaction rates -- might have a significant effect on B.

Is this concern realistic? Do I need to figure out how to have interaction between the two SIRs on an ongoing basis, or can I console myself that the non-interaction won't change things that much?

$\endgroup$
2
$\begingroup$

I can see that the combined models would still be modeling what I expect

The growth of covid-19 is not so fancy. The change in the cases are changing according to some rate $\frac{\text{d}}{\text{d}t}\log(cases)$ which is changing slowly in time.

The consequence is that you can easily fit a model to the data. People are fitting simple logistic curves or completely empirical models without any underlying mechanistic principles, such as the criticized model from the Institute for Health Metrics and Evaluation at the University of Washington.

Effectively they are more or less all the same because $\frac{\text{d}}{\text{d}t}\log(cases)$ is only changing a bit and all those models are approximately the same. Say, you could fit a quadratic curve as an approximation to all those models and it'll be the same.


but I'm a little worried about having two models with substantial populations running in parallel but not interacting after the transfer.

Your model is in principle a mechanistic model, but it is very simplified. Your model may still "work" and fit the data. But you should wonder whether the interpretation of the estimates of the epidemiological parameters will still make sense.

How to deal with this depends on whatever you wish to do with these models.

Personally I believe that the data for covid-19 (which has many flaws due to biased collection) doesn't really allow the fitting of curves and making predictions. The models should be more used for understanding the principles, and answer questions about policies in terms of relationships (e.g. 'how is the relationship between certain actions and a decrease of cases') but not with definite quantitative answers (because there is no reliable information about epidemiological parameters).


Your model could be interesting in answering the question of how the (average) parameter $\beta$ effectively changes if only a part of the population is following measures. It will help to guide questions like "What if we close churches but keep schools open?". How do the weakest links work out, is it like resistance $R_{total} = R_1 + R_2 + ...$ or is it like a piece of rope tied in series and the weakest link determines the strength of the total?

To answer such questions you should have a realistic interaction between the different nodes. The SimInf package documentation is not so clear about it, but there seems to be something like internal and external transfer which might relate respectively to transfer between individuals within the same node and transfer between individuals between different nodes.

You could also model it manually. The use of a package may not be really needed. Here is an example in python that models a cellular SIR model with interaction between the cells according to travel/commuter information.

What I expect for your model is that the growth rate will be dominated by the high $\beta$ group and that this will leak out into the low $\beta$ group effectively making the total growth rate with the same $\beta$ but not the whole population getting sick.


Example

In the example below a spatial SIR model (it is not a networked SIR model but it will give the idea) is computed where a fraction of the people are randomly selected and they have been given a 50% lower frequency of contacts. Within the normal contact group the transmission probability remains 100% of the original transmission probability, within the reduced contact group the transmission is 25% of the original, between individuals of different groups the transmission is 50% of the original.

The epidemic unfolds like an inkblob spreading out due to community transmission (which we gave a $R=2$). In addition there is a transmission to the entire plane (which we gave a smaller probability $R=0.03$), which makes that new inkblobs arise in other places.

(interesting side note: in this spatial SIR model, and also in networked SIR models, you see already early a deviation of the exponential growth, the idea of a single reproduction $R$ and basing it on the exponential growth is flawed, transmission occurs at multiple levels of depth/distance)

example of spatial SIR model

When, after 50 infections, we turn on the effect of reduced contact frequency, then you get the 'flattening of the curve'. (not just a flatter curve, but also fewer infections in total)

example of spatial SIR model with quarantine after 50 cases

The effect will be different depending on the size of the group that is following the lower contact regime (not so surprising).

effects of quarantine

Interestingly, the reduction starts of linear with the percentage of people that are following the social distancing regime. Ie. x% people that follow the regime relates to x% fewer cases of infections. But at a certain point, the drop in number of cases starts to be quick. This probably happens because the effective/average $R_0$ gets close to 1.

effects of different fractions after x days

So such a model gives an interesting insight into the mechanics of the spread. Obviously this remains a toy model that is only useful for understanding the mechanisms and understanding how policy measures may have an effect (e.g. to understand that there are non-linear and less intuitive effects, and how these effects will look like).

To get more realistic (exact) quantitative output it should be updated with a more realistic spreading across a network. That requires good information and lots of computation power. It will also still depend on a lot of guesses about epidemiological parameters. Such realism might be overkill if it the quantitative numbers remain a guesstimate. However, I believe that it remains interesting to keep some sort of interaction between the A and B group.

######
#####
#
# Spatial SIR model 
# Version 2
#
# here we make a part of the population less transmittable
#
######
######

library(progress)  # for drawing the progress bar
library("profvis") # used for optimizing the functions 

set.seed(2)

### parameters
R0 = 2     #local distribution
R1 = 0.1   #long scale distribution
L = 2*10^2

### set a fraction of people to less transmittable 
reduced = 0.5

spatialSIR2 <- function(R0 = 2, R1 = 0.1, L = 1*10^2, 
                       reduced = 0, reduction = 0.5, quarantaine_cases = 50) {

  # create LxL people in matrix
  people <- matrix(rep(0,(L)^2),L)
  # 0 indicates NOT sick


  # the algorithm will make people sick with a certain incubation time 
  # gamma distributed with mean time 4 (days)
  incubation_time <- matrix(rgamma(L^2, shape = 20, scale =4/20),L)

  # transmitancy
  # select some fraction 'reduced' of people that have reduced transferrence
  lot_trans  <-  sample(1:L^2,reduced*L^2, replace = FALSE)
  transmittance <- rep(1,L^2)
  transmittance[lot_trans] <- reduction

  # trackers for the locations of the people that got sick:
  # we start with index patient in the middle
  spread <- 0
  orderx <- L/2
  ordery <- L/2
  generation <- 1
  people[ordery,orderx] = incubation_time[ordery,orderx]

  #sick and healthy people
  # 0 = susceptible
  # 1 = sick
  # 2 = infections have spread
  sickhealthy_set <- rep(0,L^2)
  # set the index case 
  sickhealthy_set[(orderx-1)*L+ordery] = 1

  ##### details how to run the virus ######

  # compute probability density function 
  # for probabilty of spreading out to nearby locations
  Lr <- 7
  Lspread <- 1+Lr*2
  # local targets will be in a cube of LrxLr around the patient
  targets <- matrix(1:Lspread^2,Lspread)
  yt <- matrix(rep(c(1:Lspread)-(Lspread+1)/2,Lspread),Lspread)
  xt <- t(yt)
  # ps is some probability to get infected as function of distance
  ps <- c(exp(-c(Lr:1)*0.2),1,exp(-c(1:Lr)*0.2)) 
  # probs is the 2D version of ps
  probs  <- ps[xt+(Lspread+1)/2]*ps[yt+(Lspread+1)/2]  
  probs[Lspread*Lr+Lr+1] <- 0

  ### plot for visualization of the spread
  ### we uncomment this to increase spead

  #plot(orderx,ordery,xlim=c(1,L),ylim=c(1,L), 
  #     xlab = "", ylab= "",
  #     col=1,bg = 1,cex=0.2,pch=21)

  ##### run the virus ######
  # itterate all the patients in the sick_set untill all have been dealt with 
  # during this loop the number of patients increases

  sick_set <- which(sickhealthy_set == 1)

#profvis({      
#  pb <- progress_bar$new(total = L^2)
  while (0 < length(sick_set)) {
 #   pb$tick()  
    spread <- spread + 1

    # select the next first person to be sick and spread
    sick_target <- sick_set[which.min(people[sick_set])]
    sick_time <- people[sick_target]
    # coordinate of this sick person
    x <- floor((sick_target-1)/L)+1
    y <- ((sick_target-1) %% L) + 1
    g <- generation[spread]

    # selecting Rn people in the neighbourhood of the patient
    # Rn is sampled from a Poisson distribution with mean R0
    if (spread < quarantaine_cases) {
      Rn <- rpois(1,R0)
    } else {
      Rn <- rpois(1,R0*transmittance[sick_target])
    }
    if (Rn>0) {
      sel <- sample(targets,Rn, prob = probs)
      #xt[sel]
      #yt[sel]
      ## this loop picks out the R0 people 
      ## these are gonna become new patients if they are susceptible
      for (i in 1:Rn) {
        # the modulo is to patch left with right and top with bottom
        # xt,yt is the cooridinate relative to the current sick person
        # x,y is the coordinate of the current sik person
        # xq is the coordinate of the newly infected person
        xq <- (x+xt[sel[i]]-1)%%L+1  
        yq <- (y+yt[sel[i]]-1)%%L+1
        # if the 'target' is not sick yet then add it as new patient
        if  (people[yq,xq] == 0) {  
          if (spread < quarantaine_cases) {
            cont <- TRUE
          } else {
            cont <- (rbinom(1,1,transmittance[(xq-1)*L+yq])==1)
          }
          if (cont) {
            # set a sick time for the new patient
            people[yq,xq] <- sick_time + incubation_time[yq,xq]
            orderx <- c(orderx,xq)
            ordery <- c(ordery,yq)
            generation <- c(generation,g+1)
            # remove new patient from healthy set and add it to sick set
            sickhealthy_set[(xq-1)*L+yq] = 1
            sick_set <- c(sick_set,(xq-1)*L+yq)
          }
        }
      }
    }

    ### additionally make (on average) R1 random people from far away sick
    nfar <- rpois(1,R1)
    ifar <- 0
    while (ifar<nfar) {
      ifar = ifar +1
      xq <- sample(1:L,1)
      yq <- sample(1:L,1)
      ####3
      if  (people[yq,xq] == 0) {  
        if (spread < quarantaine_cases) {
          cont <- TRUE
        } else {
          cont <- (rbinom(1,1,transmittance[(xq-1)*L+yq])==1)
        }
        if (cont) {
          # set a sick time for the new patient
          people[yq,xq] <- sick_time + incubation_time[yq,xq]
          orderx <- c(orderx,xq)
          ordery <- c(ordery,yq)
          generation <- c(generation,g+1)
          # remove new patient from healthy set and add it to sick set
          sickhealthy_set[(xq-1)*L+yq] = 1
          sick_set <- c(sick_set,(xq-1)*L+yq)
        }
      }
    }

    # move patient to non-infectious group and remove from sick set
    sickhealthy_set[(x-1)*L+y] = 2          
    sick_set <- sick_set[-which(sick_set == (x-1)*L+y)]

  }
#}) profvis end
  return(list(people = people, orderx = orderx, ordery = ordery, generation = generation))
}



L = 200
set.seed(2) 
spatial2 <- spatialSIR2(L = L, R0 = 2, R1 = 0.03, reduced = 0.5)


layout(matrix(1:2,1))

# plot the epidemiological curve
times <- spatial2$people[order(spatial2$people)]
times <- times[which(times>0)]
h <- hist(times, breaks = seq(0,max(spatial2$people)+1,1), plot = FALSE )
col <- hsv(h$mids/max(spatial2$people)*0.7,0.7,1)  
plot(h$mids,h$counts, xlim = c(0,400), ylim = c(1,2000),
     xlab = "time", ylab = "newly infected",
     col=1,type = "l", log = "y",
     main="epidemiological curve")
points(h$mids,h$counts,
       col=col,bg = col,pch=21)

t <- -2:113
lines(t+1,2^(t/4)/8, lty = 2)
text(18,200, "exponential law", pos = 4 , srt = 85)

# plot the temporal spread in colours

# coordinates
ycoor <- matrix(rep(1:L,L),L)
xcoor <- t(ycoor)
# timing and colour
colvalue <- spatial2$people/max(spatial2$people)
color <- hsv(colvalue*0.7,0.7,1)  
color[spatial2$people == 0] = "white"

plot(xcoor,ycoor,xlim=c(1,L),ylim=c(1,L), 
     xlab = "", ylab= "",
     col=color,bg = color,cex=0.1,pch=21,
     main = "spatial spread of virus in time")



## computing different curves
set.seed(2)
times <- list()
pb <- progress_bar$new(total = 100)
for (i in 1:100) {
  pb$tick()
  spatial2 <- spatialSIR2(L = L, R0 = 2, R1 = 0.03, 
                          reduced = ((i-1)/20) %% 1, reduction = 0.5)
  times[[i]] <- spatial2$people[order(spatial2$people)]
  times[[i]] <- times[[i]][times[[i]]>0]
}


## computing different curves
set.seed(2)
times2 <- list()
pb <- progress_bar$new(total = 100)
for (i in 1:100) {
  pb$tick()
  spatial2 <- spatialSIR2(L = L, R0 = 2, R1 = 0.03, 
                          reduced = ((i-1)/20) %% 1, reduction = 0.75)
  times2[[i]] <- spatial2$people[order(spatial2$people)]
  times2[[i]] <- times2[[i]][times2[[i]]>0]
}


### plotting the stuff 1

plot(-100,1, xlim = c(0,500), ylim = c(1,L^2),
     xlab = "time", ylab = "cumulative infected",
     col=1,type = "l", log = "",
     main="epidemiological curves \n different fractions people \n with 50% reduced contact")

for (i in 1:100) {
  lines(times[[i]],1:length(times[[i]]), col = hsv((i%%20)/30,1,1,0.5))
  if (i %in% c(1,22,63,4,5,26,47,28,69,30,31)) {
    text(times[[i]][length(times[[i]])],length(times[[i]]), paste0(100* (((i-1)/20) %%1), " %"), 
         col = hsv((i%%20)/30,1,1), pos = 4 , cex = 0.7)
  }
}

plot(-100,1, xlim = c(0,500), ylim = c(1,6*L),
     xlab = "time", ylab = "newly infected",
     col=1,type = "l", log = "",
     main="epidemiological curves \n different fractions people \n with 50% reduced contact")

for (i in 1:100) {
  h <- hist(times[[i]], breaks = seq(0,max(times[[i]])+1,1), plot = FALSE )
  lines(h$mids,h$counts, col = hsv((i%%20)/30,1,1,0.5))
}


### plotting the stuff 2

plot(-100,1, xlim = c(0,100), ylim = c(1,L^2),
     xlab = "fraction of people with reduced contact", ylab = "cumulative infected",
     col=1,type = "l", log = "",
     main="number of infected people after x days \n Effect for different fractions people \n with 50% lower contact frequency")
max <- 0
for (i in 1:100) {
  numb1 <- sum(times[[i]]<=90)
  numb2 <- sum(times[[i]]<=120)
  numb3 <- sum(times[[i]]<=365)
  if (numb1>40) { ## not alway does the infection break out from the index case
    points(100* (((i-1)/20) %%1),numb1, 
           pch = 21, col = hsv(0,0,0), bg = hsv(0,0,0), cex = 0.7)
    points(100* (((i-1)/20) %%1),numb2,
           pch = 21, col = hsv(0,0,0), bg = hsv(0,0,0.5), cex = 0.7)
    points(100* (((i-1)/20) %%1),numb3,
           pch = 21, col = hsv(0,0,0), bg = hsv(0,0,1), cex = 0.7)
  }
  if (numb3 > max) {max <- numb3}
}

lines(c(0,100),c(max,0))
legend(0,40000, rev(c("after 90 days","after 120 days","after 365 days")),
       pch = 21, col = 1, pt.bg = rev(c(hsv(0,0,0),hsv(0,0,0.5),hsv(0,0,1))), cex = 0.7)


plot(-100,1, xlim = c(0,100), ylim = c(1,L^2),
     xlab = "fraction of people with reduced contact", ylab = "cumulative infected",
     col=1,type = "l", log = "",
     main="number of infected people after x days \n Effect for different fractions people \n with 25% lower contact frequency")
max <- 0
for (i in 1:100) {
  numb1 <- sum(times2[[i]]<=90)
  numb2 <- sum(times2[[i]]<=120)
  numb3 <- sum(times2[[i]]<=365)
  if (numb1>40) { ## not alway does the infection break out from the index case
    points(100* (((i-1)/20) %%1),numb1, 
           pch = 21, col = hsv(0,0,0), bg = hsv(0,0,0), cex = 0.7)
    points(100* (((i-1)/20) %%1),numb2,
           pch = 21, col = hsv(0,0,0), bg = hsv(0,0,0.5), cex = 0.7)
    points(100* (((i-1)/20) %%1),numb3,
           pch = 21, col = hsv(0,0,0), bg = hsv(0,0,1), cex = 0.7)
  }
  if (numb3 > max) {max <- numb3}
}

lines(c(0,200),c(max,0))
$\endgroup$
1
  • $\begingroup$ @Wayne , I have edited the answer and added a model that includes interaction between two groups with different social distancing regimes. I made the probability of transmission dependent on the product of two frequencies $f_{infected}\times f_{susceptible}$. The frequency relates to how often somebody will be showing themselves in public and the product relates to how often two individuals, one sick one susceptible, will meet. Then after the 50-th sick a fraction of the people will have their $f$ reduced by some amount (I used 50% and 25% and computed it for different fractions of people). $\endgroup$ Apr 20 '20 at 21:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.