first commit

README.md

@@ -0,0 +1,29 @@
+# Stochastic transmission model for COVID-19 transmission in a university settings
+
+## About the project
+This code accompanies the paper "High COVID-19 transmission potential associated with re-opening universities can be mitigated with layered interventions" by Ellen Brooks-Pollock, Hannah Christensen, Adam Trickey, Gibran Hemani, Emily Nixon, Amy Thomas, Katy Turner,  Adam Finn, Matt Hickman, Caroline Relton, Leon Danon, available at https://doi.org/10.1101/2020.09.10.20189696. 
+
+
+## Running the code 
+The code is written in R version 4.0.2 and C. Steps for 
+
+* Clone the repository 
+```ShellSession
+[user@host ~]$ git clone https://github.com/ellen-is/unimodel.git
+```
+* Open "UniModel.R" in R or Rstudio. You will need to install the R libraries Rcpp, tidyverse and gplots. 
+* Run script. 
+
+## Demo run and output
+The code runs ten iterations of the model with default parameters and plots the output in around 2 seconds. Increase the number of simulations with the parameter nsim, found in UniModel.R. 
+```ShellSession
+nsim=100
+```
+
+The left panel shows the number of asymptomatic and symptomatic cases over time, with one line per trajectory. The right panel shows the average (of nsim runs) number of symptomatic cases per year group. 
+
+## License
+GPL3
+
+## Contact
+[email protected]

UniModel.R

@@ -0,0 +1,81 @@
+library(Rcpp)
+library(tidyverse)
+library(gplots)
+source("plotfunctions.R") 
+
+starttime=Sys.time()
+
+
+params <- list()
+params <- within(params, {
+
+  SIMTIME=300
+  iperiod=5.0; 
+  load("data/betamat.Rdata")
+  load("data/Npop.Rdata")
+  years=read.csv("data/yearnumbers1.csv")
+  nages=length(Npop)
+
+  gam_p=1/2;
+  beta=betamat/(iperiod+1/gam_p)
+  gamma=1/iperiod;  sigma=1/(3.2);  gam_a=1/(iperiod+1/gam_p); gam_h=1/3; gam_q=1/14 ##recovery rates
+  testrate=1/2; testrate_a=rep(0,nages); testrate_a2=0; testrate_a3=0; ##testing parameters
+  eps=0.3; eps_q=0.5 ##
+  h=0.002 #hospitalisation rate
+  f=0.75 #fraction asymptomatic
+  mrate=0.038 #mortality rate
+  backgroundrate = 1e-4 #background infection rate
+
+  ## initial conditions, list within list
+  init <- within(list(), {
+    S=Npop
+    E=rep(0,nages)
+    E[sample(1:nages,10)]=1
+    A=rep(0,nages)
+    A[sample(1:nages,10)]=1
+    I=rep(0,nages)
+    I[sample(1:nages,10)]=1
+    H=rep(0,nages)
+    R=rep(0,nages)
+    R[sample(1:nages,10)]=1
+    P=rep(0,nages)
+    D=rep(0,nages)
+    Hinc=rep(0,nages)
+    Q=rep(0,nages)
+    QA=rep(0,nages)
+    N=S+E+A+I+H+R+P+Q+QA
+  })
+})
+attach(params)
+set.seed(params$seed)
+sourceCpp('studentmodel_reactive2.cpp')
+
+nsim=10
+par(mar=c(5,5,1,1),mfrow=c(1,2))
+
+result.rep=c()
+for(nn in 1:nsim)
+{
+  set.seed(nn)
+  result <- c19uni(params)
+  result$nsim <- nn
+  result$time[params$SIMTIME]=(params$SIMTIME-1)
+  result.rep=rbind(result.rep,result)
+}
+
+
+Infecteds=getinfecteds(result.rep,years)
+allinf=getasymp(result.rep,years)
+
+Infecteds %>%
+  group_by(time) %>%
+  summarise_all(mean) -> meanout
+
+maxbaseline=max(allinf$total)*1.2
+plotcases(allinf,Infecteds,maxbaseline,'baseline')
+plotyears(meanout)
+
+
+end.time = Sys.time()
+
+print(end.time-starttime)

data/yearnumbers1.csv

@@ -0,0 +1,154 @@
+"x"
+"0"
+"0"
+"0"
+"0"
+"1"
+"1"
+"1"
+"1"
+"1"
+"1"
+"1"
+"1"
+"1"
+"1"
+"1"
+"1"
+"1"
+"1"
+"1"
+"1"
+"1"
+"1"
+"1"
+"1"
+"1"
+"1"
+"1"
+"1"
+"1"
+"1"
+"1"
+"2"
+"2"
+"2"
+"2"
+"2"
+"2"
+"2"
+"2"
+"2"
+"2"
+"2"
+"2"
+"2"
+"2"
+"2"
+"2"
+"2"
+"2"
+"2"
+"2"
+"2"
+"2"
+"2"
+"2"
+"2"
+"2"
+"3"
+"3"
+"3"
+"3"
+"3"
+"3"
+"3"
+"3"
+"3"
+"3"
+"3"
+"3"
+"3"
+"3"
+"3"
+"3"
+"3"
+"3"
+"3"
+"3"
+"3"
+"3"
+"3"
+"3"
+"3"
+"3"
+"4"
+"4"
+"4"
+"4"
+"4"
+"4"
+"4"
+"4"
+"4"
+"4"
+"4"
+"4"
+"4"
+"4"
+"4"
+"4"
+"4"
+"4"
+"4"
+"4"
+"4"
+"5"
+"5"
+"5"
+"5"
+"R"
+"R"
+"R"
+"R"
+"R"
+"R"
+"R"
+"R"
+"R"
+"R"
+"R"
+"R"
+"R"
+"R"
+"R"
+"R"
+"R"
+"R"
+"R"
+"R"
+"R"
+"R"
+"R"
+"T"
+"T"
+"T"
+"T"
+"T"
+"T"
+"T"
+"T"
+"T"
+"T"
+"T"
+"T"
+"T"
+"T"
+"T"
+"T"
+"T"
+"T"
+"T"
+"T"
+"T"
+"T"

plotfunctions.R

@@ -0,0 +1,86 @@
+
+getinfecteds = function(result.rep,years)
+{
+  Infecteds = result.rep[,c((3*nages+2):((3+1)*nages+1),1,length(result.rep[1,]))]
+  Infecteds$year1 = rowSums(Infecteds[which(years=="1")])
+  Infecteds$year2 = rowSums(Infecteds[which(years=="2")])
+  Infecteds$year3 = rowSums(Infecteds[which(years=="3")])
+  Infecteds$year4 = rowSums(Infecteds[which(years=="4")])
+  Infecteds$yearR = rowSums(Infecteds[which(years=="R")])
+  Infecteds$yearT = rowSums(Infecteds[which(years=="T")])
+  Infecteds$total = rowSums(Infecteds[1:nages])
+  return(Infecteds)
+}
+getasymp = function(result.rep,years)
+{
+
+  allinf = result.rep[,c((2*nages+2):((2+1)*nages+1),
+                         (6*nages+2):((6+1)*nages+1),1,length(result.rep[1,]))]
+
+  allinf$total = rowSums(allinf[1:(length(allinf[1,])-2)])
+
+  return(allinf)
+}
+getallinf = function(result.rep,years)
+{
+
+  allinf = result.rep[,c((1*nages+2):((1+1)*nages+1),
+                         (2*nages+2):((2+1)*nages+1),
+                         (3*nages+2):((3+1)*nages+1),
+                         (6*nages+2):((6+1)*nages+1),
+                         (8*nages+2):((8+1)*nages+1),
+                         (9*nages+2):((9+1)*nages+1),
+                         1,length(result.rep[1,]))]
+  allinf$total = rowSums(allinf[1:(length(allinf[1,])-2)])
+
+  return(allinf)
+}
+
+
+plotyears = function(meanout)
+{
+  studentcolsa=hcl.colors(10, "YlOrRd", rev = FALSE)
+  studentcols=hcl.colors(10, "YlOrRd", rev = FALSE,alpha=0.1)
+  startofterm=as.Date("28/09/2020",format="%d/%m/%Y")
+  Xmas=as.Date("21/12/2020",format="%d/%m/%Y")-startofterm
+  days=1:length(meanout$year1)
+
+  studentcolsa=hcl.colors(10, "YlOrRd", rev = FALSE)
+  studentcols=hcl.colors(10, "YlOrRd", rev = FALSE,alpha=0.1)
+
+  plot(days+startofterm,meanout$year1,type="l",lwd=2,col=studentcolsa[1],
+       xlab="",cex.axis=1.1,cex.lab=1.3,ylab="# symptomatic")
+
+  lines(days+startofterm,meanout$year2,type="l",lwd=2,col=studentcolsa[2])
+  lines(days+startofterm,meanout$year3,type="l",lwd=2,col=studentcolsa[3])
+  lines(days+startofterm,meanout$year4,type="l",lwd=2,col=studentcolsa[4])
+  lines(days+startofterm,meanout$yearT,type="l",lwd=2,col=studentcolsa[5])
+  lines(days+startofterm,meanout$yearR,type="l",lwd=2,col=studentcolsa[6])
+  legend('topright',c(paste('year',1:4),"PGT","PGR"),col=studentcolsa,lwd=3,cex=1.0,box.col=NA)
+
+  abline(v=as.Date("21/12/2020",format="%d/%m/%Y"),lty=2)
+}
+
+
+plotcases = function(allinf,Infecteds,maxbaseline,stratname)
+{
+  studentcolsa=hcl.colors(10, "YlOrRd", rev = FALSE,alpha=0.6)
+  studentcols=hcl.colors(10, "YlOrRd", rev = FALSE,alpha=0.1)
+  startofterm=as.Date("28/09/2020",format="%d/%m/%Y")
+  Xmas=as.Date("21/12/2020",format="%d/%m/%Y")-startofterm
+  n=1
+  days=1:length(allinf$total[allinf$nsim==n])
+  plot(startofterm+days,allinf$total[allinf$nsim==n],type="l",lwd=2,col=rgb(0,0.5,0,0.1),
+       xlab="",cex.axis=1.1,cex.lab=1.3,
+       ylim=c(1,1.0*maxbaseline),ylab="# infected",log="")
+
+  for(n in 1:nsim)
+  {
+    lines(startofterm+days,Infecteds$total[Infecteds$nsim==n],type="l",lwd=2,col=studentcols[1])
+
+    lines(startofterm+days,allinf$total[allinf$nsim==n],type="l",lwd=2,col=rgb(0,0.5,0,0.1))
+  }
+  abline(v=as.Date("21/12/2020",format="%d/%m/%Y"),lty=2)
+  legend('topright',c("symptomatic","asymptomatic"),col=c(studentcolsa[1],rgb(0,0.5,0,0.5)),pch=19,box.col = NA)
+}
+