Laying out the maths

The model is a simple stratified SIR. The principle of parsomony applies. The aim is to make the model as simple as possible, but not more so.

The population is divided into a vulnerable fraction who may contract the illness with some risk of death and the fraction who contract the illness and have a very high probability of complete recover

The same model is run through in two phases. The first phase represents return to normality for those who are not vulnerable, but partial shielding for the vulnerable. The second represents complete return to normality.

‘Safe’ fraction

\(\frac {dS_1}{dt} = - \beta_1 \frac{S_1I_1}{N_2}\)

\(\frac {dI_1}{dt} = \beta_1 \frac{S_1I_1}{N_1}- \gamma \frac{I_1}{N_1}\)

\(\frac {dR_1}{dt} = \gamma \frac{I_1}{N_1}\)

‘Vulnerable’ fraction

This really needs thought!

\(\frac {dI_2}{dt} = \beta_2 \frac{S_2I_2}{N} + \beta_{12} \frac{S_2I_1}{N_2} - \gamma \frac{I_2}{N_2}\)

\(\frac {dS_2}{dt} = - \beta_2 \frac{S_2I_1}{N_2} - \beta_{12} \frac{S_2I_1}{N_2}\)

\(\frac {dR_2}{dt} = \gamma \frac{I_2}{N_1}\)

Mortality is then simply calculated by mutiplying the final number of removed vulnerables by the IFR (note the same simple calculation can be used for mortality of the less vulnerable is needed)

\(IFR_{vulnerable} R_2N_2 + IFR_{safe} R_1N_1\)

So the idea is that there is only one parameter \(\gamma\) for recovery … which could be realxed if really needed to become \(\gamma_1\) and \(\gamma_2\). But there are three parameters for transmission

\(\beta_1\) for transmission within the “safe” fraction

\(\beta_2\) for transmission within the “vulnerable” fraction and \(\beta_{12}\) for transmission between the safe and vulnerable fraction

With this element of the equation representing the transfer …

\(\beta_{12} \frac{S_2I_1}{N_2}\)

Next step … run phase one for x period of days starting with the situation estimated as the present state fo the population. Then reparamramterise the whole thing for the second period, initiaiting the compartments where we left off at the end of the first.

## Create an GBD function

gbd <- function(time, state, parameters) {
  
  with(as.list(c(state, parameters)), {
    
 
    dS1 <- -beta1 * S1 * I1
    dI1 <-  beta1 * S1 * I1 - gamma * I1
    dR1 <- gamma * I1 
    
    dS2 <- - (beta2 * S2 * I2 + beta12 * S2 * I1)
    dI2 <-  (beta2 * S2 * I2) + (beta12 * S2 * I1) - (gamma * I2)
    dR2 <- gamma * I2 
    
    
    return(list(c(dS1, dI1, dR1,dS2, dI2, dR2 )))
  })
}

run_gbd<-function(
  beta1=0.5,beta2=0.2 , gamma= 0.1,
  S1=0.9,I1=0.1, R1=0,
  S2=0.99,I2=0.01, R2=0,
  beta12 =0.01,
  days=60
  )
{
parameters <- c(beta1 = beta1,beta2=beta2, beta12=beta12,gamma=gamma)
init<- c(S1 = S1, I1 = I1,R1 = R1, S2 = S2, I2 = I2,R2 = R2)
times<- seq(0, days, by = 1)
out <- ode(y = init, times = times, func = gbd, parms = parameters)
d <- as.data.frame(out)        
d
}
sd<-run_gbd(beta2=0, beta12=0.01,beta1=0.6)
dygraph(sd) 
ds<-pivot_longer(sd,cols=2:7)


ds %>% filter(name %in% c("R1","R2")) %>%
ggplot(aes(x=time,y=value,color=name)) + geom_line()

ds %>% filter(name %in% c("I1","I2")) %>%
ggplot(aes(x=time,y=value,color=name)) + geom_line()

final_R2<-function(R01=3,R02=1,R012=0.5,gamma=0.2,days=300){
beta1<-R01*gamma
beta2<-R02*gamma
beta12<-R012*gamma
sd<-run_gbd(beta2, beta12,beta1)
max(sd$R2) * 100
}


d<-data.frame(R01=rep(1:10,each=10),R02=rep(1:10,by=10),R012=runif(100,0.2,2.2),gamma=0.2)

d %>% rowwise %>% mutate(res=final_R2(R01,R02,R012)) %>% aqm::dt()