diff --git a/README.md b/README.md
new file mode 100644
index 0000000..049dbb6
--- /dev/null
+++ b/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
+Ellen.Brooks-Pollock@bristol.ac.uk
diff --git a/UniModel.R b/UniModel.R
new file mode 100644
index 0000000..5aa2f72
--- /dev/null
+++ b/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)
diff --git a/data/Npop.Rdata b/data/Npop.Rdata
new file mode 100644
index 0000000..efeb4e6
Binary files /dev/null and b/data/Npop.Rdata differ
diff --git a/data/betamat.RData b/data/betamat.RData
new file mode 100644
index 0000000..174eb3b
Binary files /dev/null and b/data/betamat.RData differ
diff --git a/data/yearnumbers1.csv b/data/yearnumbers1.csv
new file mode 100644
index 0000000..59f2302
--- /dev/null
+++ b/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"
diff --git a/plotfunctions.R b/plotfunctions.R
new file mode 100644
index 0000000..e45f19c
--- /dev/null
+++ b/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)
+}
+
diff --git a/studentmodel_reactive2.cpp b/studentmodel_reactive2.cpp
new file mode 100755
index 0000000..a07c266
--- /dev/null
+++ b/studentmodel_reactive2.cpp
@@ -0,0 +1,440 @@
+#include <Rcpp.h>
+using namespace Rcpp;
+
+// [[Rcpp::export]]
+List c19uni(List params) {
+  
+  int j=0, J=0;
+  double k=0.0, m=0.0, i=0.0;
+  //printf("1\n");
+  // pull out list w/in list into its own object
+  List init = params["init"];
+  List S0 = init["S"];
+  List E0 = init["E"];
+  List A0 = init["A"];
+  List P0 = init["P"];
+  List I0 = init["I"];
+  List R0 = init["R"];
+  List H0 = init["H"];
+  List Q0 = init["Q"];
+  List QA0 = init["QA"];
+  List D0 = init["D"];
+  List N0 = init["N"];
+  
+  //printf("2\n");
+  // use Rcpp as() function to "cast" R vector to cpp scalar
+  int nsteps = as<int>(params["SIMTIME"]);
+  int nages = as<int>(params["nages"]);
+  
+  //printf("3\n");
+  // initialize each state vector in its own vector
+  NumericMatrix SS(nsteps,nages);
+  NumericMatrix EE(nsteps,nages);
+  NumericMatrix AA(nsteps,nages);
+  NumericMatrix PP(nsteps,nages);
+  NumericMatrix II(nsteps,nages);
+  NumericMatrix HH(nsteps,nages);
+  NumericMatrix RR(nsteps,nages);
+  NumericMatrix QQ(nsteps,nages);
+  NumericMatrix QA(nsteps,nages);
+  NumericMatrix DD(nsteps,nages);
+  NumericMatrix NN(nsteps,nages);
+
+  //printf("4\n");
+  for(j=0; j<nages; j++)
+ {
+     SS(0,j) = S0[j];
+    EE(0,j) = E0[j];
+  AA(0,j) = A0[j];
+ PP(0,j) = P0[j];
+  II(0,j) = I0[j];
+  HH(0,j) = H0[j];
+ RR(0,j) = R0[j];
+ QQ(0,j) = Q0[j];
+ QA(0,j) = QA0[j];
+ DD(0,j) = D0[j];
+ NN(0,j) = N0[j];
+  }
+  //printf("5\n");
+  // fill time w/zeros
+  NumericVector time(nsteps);
+  //printf("5B\n");
+  // pull out params for easy reading
+  //[1] "init"       "mrate"      "f"          "h"          "eps"        "testrate_a" "testrate"   "gam_q"      "gam_p"      "gam_h"      "gam_a"     
+  //[12] "sigma"      "gamma"      "beta1"      "rzero"      "nages"      "years"      "groupnames" "Npop"       "iperiod"    "SIMTIME"    "seed"
+  double sigma = params["sigma"];
+  double mrate = params["mrate"];
+  double f1 = params["f"];
+  double eps = params["eps"];
+  double eps_q = params["eps_q"];
+  double h = params["h"];
+  double gamma = params["gamma"];
+  double gam_p = params["gam_p"];
+  double gam_h = params["gam_h"];
+  double gam_a = params["gam_a"];
+  double gam_q = params["gam_q"];
+  double backgroundrate = params["backgroundrate"];
+  
+  double testrate = params["testrate"];
+  double testrate_a2 = params["testrate_a2"];
+  double testrate_a3 = params["testrate_a3"];
+  //printf("5C\n");
+  NumericVector testrate_a = params["testrate_a"];
+  NumericMatrix beta = params["beta"];
+  NumericVector testrate_a_use = testrate_a;
+  
+  //printf("6\n");
+  // declare transition vectors
+  NumericVector SE(nages); //set later
+  NumericVector EA(nages,f1 * sigma);
+  NumericVector EP(nages,(1.0 - f1) * sigma);
+  NumericVector PInf(nages,gam_p);
+  NumericVector AR(nages,gam_a);
+  NumericVector IH(nages,h * gamma);
+  NumericVector IR(nages,(1.0 - h) * gamma);
+  NumericVector HR(nages,(1 - mrate) * gam_h);
+  NumericVector HD(nages,mrate * gam_h);
+  NumericVector IQ(nages,testrate);
+  NumericVector AQ(nages);//set later
+  //NumericVector AQ testrate_a;//set later
+  NumericVector PQ(nages);//set later
+  NumericVector QR(nages,gam_q);
+  
+  double prob = 0.0;
+
+  int flag=0,daycount=1,newcases=0,newcases_then=0,delta=0; 
+  double Etot=0.0,Ntot=0.0,Ntot2=0.0;
+  for (int t= 0; t < (nsteps-1); t++) {
+    time[t] = t;
+    
+    for(j = 0; j < nages; j++)
+    {
+      SS(t+1,j) = SS(t,j);
+      EE(t+1,j) = EE(t,j);
+      AA(t+1,j) = AA(t,j);
+      PP(t+1,j) = PP(t,j);
+      II(t+1,j) = II(t,j);
+      RR(t+1,j) = RR(t,j);
+      QQ(t+1,j) = QQ(t,j);
+      QA(t+1,j) = QA(t,j);
+      HH(t+1,j) = HH(t,j);
+      DD(t+1,j) = DD(t,j);
+      
+      if(EE(t+1,j)<0){printf("set up E[%i,%i]=%f\n",t+1,j,EE(t+1,j)); break;}
+    }
+    Etot=0;
+    for(j = 0; j < nages; j++) Etot = Etot + EE(t+1,j);
+    //printf("t=%i Etot=%f\n",t+1,Etot);
+    if(isnan(Etot)){break;}
+    
+    Ntot=0.0; Ntot2=0.0; 
+    for(j = 0; j < nages; j++) 
+    {
+      Ntot2 = SS(t+1,j) + EE(t+1,j) + AA(t+1,j) + PP(t+1,j) + II(t+1,j) + RR(t+1,j) + HH(t+1,j) + QQ(t+1,j) + QA(t,j);
+      Ntot = SS(t,j) + EE(t,j) + AA(t,j) + PP(t,j) + II(t,j) + RR(t,j) + HH(t,j) + QQ(t,j) + QA(t,j);
+    }
+    if(Ntot!=Ntot2){printf("1 t=%i N1 %f N2 %f\n",t, Ntot,Ntot2); break;}
+    
+    flag=0; 
+    //printf("t = %i\n",t);
+    for(J = 0; J < nages; J++) 
+      {
+        delta=0;
+        SE[J] = 0.0;
+        prob=0.0;
+        for(j = 0; j < nages; j++)
+        {
+          if(j==J){delta=1;}
+          if(NN(t,j)>0)
+          {
+            //prob += beta(J, j) * (eps_q*eps*QA(t,j) + eps_q*QQ(t,j) + eps*AA(t,j)+PP(t,j)+II(t,j)) / NN(t,j);
+            if(eps_q==0)
+            {
+              if(t>=70 & t<=84)
+              {
+              prob += beta(J, j) * (eps_q*delta*eps*QA(t,j) + eps_q*delta*QQ(t,j) + eps*AA(t,j)+PP(t,j)+II(t,j)) / NN(t,j);
+              }
+              else{
+                prob += beta(J, j) * (eps_q*delta*eps*QA(t,j) + eps_q*delta*QQ(t,j) + delta*eps*AA(t,j)+PP(t,j)+delta*II(t,j)) / NN(t,j);
+              }
+            }
+            else{
+              prob += beta(J, j) * (eps_q*delta*eps*QA(t,j) + eps_q*delta*QQ(t,j) + eps*AA(t,j)+PP(t,j)+II(t,j)) / NN(t,j);
+            }
+            
+          }
+          
+          //printf("J %i j %i beta %f prob %f NN(t,j) %f\n",J,j,beta(J,j),prob,NN(t,j));
+        }
+        
+        SE[J] = 1 - exp(-prob - backgroundrate);
+    }
+
+    // set AQ
+    for(j = 0; j < nages; j++) AQ[j] = testrate_a_use[j];
+    // set PQ
+    for(j = 0; j < nages; j++) PQ[j] = testrate_a_use[j];
+
+    
+    /////////////////////////
+    // State Equations
+    /////////////////////////
+    // discrete-time model
+    // SE
+    for(j = 0; j < nages; j++) 
+      {
+        i=0.0;
+        prob = SE[j];
+        //printf("num susceptibles S[%i,%i]=%f prob %f SE %f\n",t,j,SS(t,j),prob,SE[j]);
+        i = R::rbinom(SS(t,j), prob);
+        //printf("new infections i=%f\n",i);
+        SS(t+1,j) = SS(t+1,j) - i;
+        //EE(t+1,j) = EE(t+1,j) + i;
+        EE(t+1,j) += i;
+
+      } 
+    
+    Ntot=0.0; Ntot2=0.0; 
+    for(j = 0; j < nages; j++) 
+    {
+      Ntot2 = SS(t+1,j) + EE(t+1,j) + AA(t+1,j) + PP(t+1,j) + II(t+1,j) + RR(t+1,j) + HH(t+1,j) + QQ(t+1,j) + QA(t+1,j);
+      Ntot = SS(t,j) + EE(t,j) + AA(t,j) + PP(t,j) + II(t,j) + RR(t,j) + HH(t,j) + QQ(t,j) + QA(t,j);
+    }
+    if(Ntot!=Ntot2){printf("2 t=%i N1 %f N2 %f\n",t, Ntot,Ntot2); break;}
+    
+    // EA or EP
+    for(j = 0; j < nages; j++) {
+      i=0.0; k=0;
+      if(EE(t,j)>0)
+      {
+        prob = 1 - exp(-(EA[j] - EP[j]));
+        i = R::rbinom(EE(t,j), prob);
+        if(i > 0) {
+          prob = EA[j] / (EA[j] + EP[j]);
+          k = R::rbinom(i, prob);
+        }
+      }
+      if(isnan(EE(t+1,j))){break;}
+      //EE(t+1,j) = EE(t+1,j) - i;
+      EE(t+1,j) -= i;
+      AA(t+1,j) = AA(t+1,j) + k;
+      PP(t+1,j) = PP(t+1,j) + i - k;
+      if(EE(t+1,j)<0){printf("i= %f E[%i,%i]=%f, E[%i,%i]=%f\n",i,t,j,EE(t,j),t+1,j,EE(t+1,j)); break;}
+      //printf("2. E[%i,%i]=%f\n",t+1,j,EE(t+1,j));
+    }
+    
+    
+    Ntot=0.0; Ntot2=0.0; 
+    for(j = 0; j < nages; j++) 
+    {
+      Ntot2 = SS(t+1,j) + EE(t+1,j) + AA(t+1,j) + PP(t+1,j) + II(t+1,j) + RR(t+1,j) + HH(t+1,j) + QQ(t+1,j) + QA(t+1,j);
+      Ntot = SS(t,j) + EE(t,j) + AA(t,j) + PP(t,j) + II(t,j) + RR(t,j) + HH(t,j) + QQ(t,j) + QA(t,j);
+    }
+    if(Ntot!=Ntot2){printf("3 t=%i N1 %f N2 %f\n",t, Ntot,Ntot2); break;}
+    
+    
+    // AR or A-QA
+    for(j = 0; j < nages; j++) {
+      i=0; k=0.0;
+      if(AA(t,j)>0)
+      {
+        prob = 1 - exp(-AR[j] - AQ[j]);
+        i = R::rbinom(AA(t,j), prob);
+        if(i>0)
+        {
+          prob = AR[j] / (AR[j]  + AQ[j]);
+          k = R::rbinom(i, prob);
+          //printf("k=%f\n",k);
+        }
+        AA(t+1,j) = AA(t+1,j) - i;
+        RR(t+1,j) = RR(t+1,j) + k;
+        //QQ(t+1,j) = QQ(t+1,j) + i - k;
+        QA(t+1,j) = QA(t+1,j) + i - k;
+        newcases = newcases +i -k;
+      }
+
+    }
+    
+    
+    Ntot=0.0; Ntot2=0.0; 
+    for(j = 0; j < nages; j++) 
+    {
+      Ntot2 = SS(t+1,j) + EE(t+1,j) + AA(t+1,j) + PP(t+1,j) + II(t+1,j) + RR(t+1,j) + HH(t+1,j) + QQ(t+1,j) + QA(t+1,j);
+      Ntot = SS(t,j) + EE(t,j) + AA(t,j) + PP(t,j) + II(t,j) + RR(t,j) + HH(t,j) + QQ(t,j) + QA(t,j);
+    }
+    if(Ntot!=Ntot2){printf("4 t=%i N1 %f N2 %f\n",t, Ntot,Ntot2); break;}
+    
+    
+    
+    // PI or PQ
+    for(j = 0; j < nages; j++) {
+      i=0; k=0;
+      prob = 1 - exp(-PInf[j] - PQ[j]);
+      i = R::rbinom(PP(t,j), prob);
+      if(i>0)
+      {
+        prob = PInf[j] / (PInf[j]  + PQ[j]);
+        k = R::rbinom(i, prob);
+      }
+      PP(t+1,j) = PP(t+1,j) - i;
+      II(t+1,j) = II(t+1,j) + k;
+      QQ(t+1,j) = QQ(t+1,j) + i - k;
+      newcases = newcases + i - k;
+    }
+    
+    
+    Ntot=0.0; Ntot2=0.0; 
+    for(j = 0; j < nages; j++) 
+    {
+      Ntot2 = SS(t+1,j) + EE(t+1,j) + AA(t+1,j) + PP(t+1,j) + II(t+1,j) + RR(t+1,j) + HH(t+1,j) + QQ(t+1,j) + QA(t+1,j);
+      Ntot = SS(t,j) + EE(t,j) + AA(t,j) + PP(t,j) + II(t,j) + RR(t,j) + HH(t,j) + QQ(t,j) + QA(t,j);
+    }
+    if(Ntot!=Ntot2){printf("5 t=%i N1 %f N2 %f\n",t, Ntot,Ntot2); break;}
+    
+    
+    // IH or IR or IQ
+    for(j = 0; j < nages; j++) {
+      i=0; k=0; m=0;
+      prob = 1 - exp(-(IH[j] + IR[j] + IQ[j]));
+      i = R::rbinom(II(t,j), prob);
+      if(i > 0) {
+        prob = IH[j] / (IH[j] + IR[j] + IQ[j]);
+        k = R::rbinom(i, prob);
+        if((i-k) > 0)
+        {
+          prob = IR[j] / (IR[j] + IQ[j]);
+          m = R::rbinom(i-k, prob);
+          //printf("m=%f\n",m);
+        }
+      }
+      II(t+1,j) = II(t+1,j) - i;
+      HH(t+1,j) = HH(t+1,j) + k;
+      RR(t+1,j) = RR(t+1,j) + m;
+      QQ(t+1,j) = QQ(t+1,j) + i - k - m;
+      newcases = newcases + i -k -m;
+    }
+
+    
+    Ntot=0.0; Ntot2=0.0; 
+    for(j = 0; j < nages; j++) 
+    {
+      Ntot2 = SS(t+1,j) + EE(t+1,j) + AA(t+1,j) + PP(t+1,j) + II(t+1,j) + RR(t+1,j) + HH(t+1,j) + QQ(t+1,j) + QA(t+1,j);
+      Ntot = SS(t,j) + EE(t,j) + AA(t,j) + PP(t,j) + II(t,j) + RR(t,j) + HH(t,j) + QQ(t,j) + QA(t,j);
+    }
+    if(Ntot!=Ntot2){printf("6 t=%i N1 %f N2 %f\n",t, Ntot,Ntot2); break;}
+    
+    // HR or HD
+    for(j = 0; j < nages; j++) {
+      i=0; k=0;
+      prob = 1 - exp(-(HR[j] + HD[j]));
+      i = R::rbinom(HH(t,j), prob);
+      if(i > 0) {
+        prob = HR[j] / (HR[j] + HD[j]);
+        k = R::rbinom(i, prob);
+      }
+      HH(t+1,j) = HH(t+1,j) - i;
+      RR(t+1,j) = RR(t+1,j) + i;
+      DD(t+1,j) = DD(t+1,j) +i -k;
+    }
+    
+    Ntot=0.0; Ntot2=0.0; 
+    for(j = 0; j < nages; j++) 
+    {
+      Ntot2 = SS(t+1,j) + EE(t+1,j) + AA(t+1,j) + PP(t+1,j) + II(t+1,j) + RR(t+1,j) + HH(t+1,j) + QQ(t+1,j) + QA(t+1,j);
+      Ntot = SS(t,j) + EE(t,j) + AA(t,j) + PP(t,j) + II(t,j) + RR(t,j) + HH(t,j) + QQ(t,j) + QA(t,j);
+    }
+    if(Ntot!=Ntot2){printf("7 t=%i N1 %f N2 %f\n",t, Ntot,Ntot2); break;}
+    
+    // QR
+    for(j = 0; j < nages; j++) {
+      if(QQ(t,j)>0)
+      {
+        prob = 1 - exp(-QR[j]);
+        i = R::rbinom(QQ(t,j), prob);
+        QQ(t+1,j) = QQ(t+1,j) - i;
+        RR(t+1,j) = RR(t+1,j) + i;
+      }
+    }
+
+    for(j = 0; j < nages; j++) {
+      if(QA(t,j)>0)
+      {
+        prob = 1 - exp(-QR[j]);
+        i = R::rbinom(QA(t,j), prob);
+        QA(t+1,j) = QA(t+1,j) - i;
+        RR(t+1,j) = RR(t+1,j) + i;
+      }
+    }
+    
+    
+    Ntot=0.0; Ntot2=0.0; 
+    for(j = 0; j < nages; j++) 
+    {
+      Ntot2 = SS(t+1,j) + EE(t+1,j) + AA(t+1,j) + PP(t+1,j) + II(t+1,j) + RR(t+1,j) + HH(t+1,j) + QQ(t+1,j) + QA(t+1,j);
+      Ntot = SS(t,j) + EE(t,j) + AA(t,j) + PP(t,j) + II(t,j) + RR(t,j) + HH(t,j) + QQ(t,j) + QA(t,j);
+    }
+    if(Ntot!=Ntot2){printf("8 t=%i N1 %f N2 %f\n",t, Ntot,Ntot2); break;}
+    
+    for(j = 0; j < nages; j++) 
+    {
+      NN(t+1,j) = SS(t+1,j) + EE(t+1,j) + AA(t+1,j) + PP(t+1,j) + II(t+1,j) + RR(t+1,j) + HH(t+1,j) + QQ(t+1,j) + QA(t+1,j);
+      
+      if(isnan(NN(t+1,j)) || EE(t+1,j)<0)
+      {
+      printf("S=%f E=%f A=%f P=%f I=%f R=%f Q=%f H=%f N=%f\n",SS(t+1,j),EE(t+1,j),AA(t+1,j),PP(t+1,j),II(t+1,j),RR(t+1,j),QQ(t+1,j),HH(t+1,j),NN(t+1,j));
+      break;
+      }
+    }
+    
+    
+    // for(j = 0; j < nages; j++){
+    //   if(II(t+1,j)>(II(t,j)+5)){
+    //     testrate_a_use[j] = testrate_a2;
+    //     printf("t=%i j=%i II[t+1,j]=%f II[t,j]=%f, increase testing\n",t+1,j,II(t+1,j),II(t,j));
+    //     }
+    //   else{testrate_a_use[j] = testrate_a3;}
+    // } 
+    
+    daycount++;
+    if(daycount==7)
+    {
+      if(newcases>0)
+      {
+        if((newcases/(newcases_then+1)) >= 1)
+        {
+            for(j = 0; j < nages; j++){
+              testrate_a_use[j] = testrate_a2;
+            }
+            //printf("t=%i newcases %i newcasesold %i, increase testing\n",t+1,newcases,newcases_then);
+        }
+        else{
+          for(j = 0; j < nages; j++){
+            testrate_a_use[j] = testrate_a3;
+          }
+          //printf("t=%i newcases %i newcasesold %i, decrease testing\n",t+1,newcases,newcases_then);
+        }
+      }
+      newcases_then = newcases;
+      newcases=0;
+      daycount=1;
+    }
+    
+  };//end of the t loop
+  
+  // Return results as data.frame
+  DataFrame sim = DataFrame::create(
+    Named("time") = time,
+    Named("S") = SS,
+    Named("E") = EE,
+    Named("A") = AA,
+    Named("I") = II,
+    Named("R") = RR,
+    Named("H") = HH,
+    Named("P") = PP,
+    Named("D") = DD,
+    Named("Q") = QQ,
+    Named("QA") = QA,
+    Named("N") = NN
+  );
+  return sim;
+};