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)
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)
The effect will be different depending on the size of the group that is following the lower contact regime (not so surprising).
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.
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))
