R/shiny2strain.r

In epimdr2: Functions and Data for "Epidemics: Models and Data in R (2nd Edition)"

#' Launch a shiny-app simulating a two-strain SIR model
#' @details
#' Launch app for details
#' @examples
#' if(interactive()){twostrain.app}
#' @export

twostrain.app=shinyApp(
ui = pageWithSidebar(
  headerPanel("A two-strain SIR model"),
  sidebarPanel(
    sliderInput("beta1", "Transmission 1 (yr^-1):", 500,
                min = 0, max = 3000),
    sliderInput("beta2", "Transmission 2 (yr^-1):", 750,
                min = 0, max = 3000),
    sliderInput("Ip", "Infectious period (days)", 5,
                min = 1, max = 100),
    sliderInput("mu", "birth rate (per year):", 0.02,
                min = 0, max = .1),
    sliderInput("theta", "Theta:", 0.5,
                min = 0, max = 1),
    sliderInput("xi", "Xi:", 0.5,
                min = 0, max = 1),
    sliderInput("pi", "Pi:", 0.5,
                min = 0, max = 1),
    sliderInput("T", "Time range:",
                min = 0, max = 30, value = c(10,30)),
    #checkboxInput("lg", "Log", FALSE), width=3),
    checkboxInput("lg", "4th root", FALSE), width=3),
  mainPanel(
    tabsetPanel(
      tabPanel("Parameters", dataTableOutput("table1")), 
      tabPanel("Time", plotOutput("plot1")),
      tabPanel("S* plot", plotOutput("plot2")),
      tabPanel("Phase plane", plotOutput("plot3")),
      tabPanel("Details", 
               withMathJax(
                 helpText("MODEL:"),
                 helpText("Susceptible $$\\frac{dS}{dt} = \\mu (N - S) - \\phi$$"),
                 helpText("$$\\phi = \\frac{\\beta_1 * I1 +\\beta_2 * I2+\\Theta * (\\beta_1 * J1+\\beta_2* J_2)}{N}$$"),
                 helpText("Strain 1 I $$\\frac{dI_1}{dt} = \\frac{(\\beta_1 * I_1 + \\Theta * \\beta_1 * J_1) * S}{N} - (\\gamma + \\mu) * I_1$$"),
                 helpText("Strain 2 I $$\\frac{dI_2}{dt} = \\frac{(\\beta_2 * I_2 + \\Theta * \\beta_2 * J_2) * S}{N} - (\\gamma + \\mu) * I_2$$"),
                 helpText("Immune to 1 only $$\\frac{dR_1}{dt} = \\Pi * \\gamma * I_1 - (\\beta_2 * I_2 + \\Theta * \\beta_2 * J_2) * \\Xi * R_1 / N -  \\mu * R_1$$"),
                 helpText("Immune to 2 only $$\\frac{dR_2}{dt} = \\Pi * \\gamma * I_2 - (\\beta_1 * I_1 + \\Theta * \\beta_1 * J_1) * \\Xi * R_2 / N -  \\mu * R_2$$"),
                 helpText("Strain 1 J $$\\frac{dJ_1}{dt} = (\\beta_1 * I_1 + \\Theta * \\beta_1 * J_1) * \\Xi * R_2 / N - \\gamma * J_1 - \\mu * J_1$$"),
                 helpText("Strain 2 J $$\\frac{dJ_2}{dt} = (\\beta_2 * I_2 + \\Theta * \\beta_2 * J_2) * \\Xi * R_1 / N - \\gamma * J_2 - \\mu * J_2$$"),
                 helpText("Removed $$\\frac{dR}{dt} = (1-\\Pi) * \\gamma * (I_1 + I_2) + \\gamma * (J_1 + J_2) - \\mu * R$$"),  
                 helpText("REFERENCE: ")
               ))
      
    )
  )
),


#R_{0,1}=&\frac{\beta}{\gamma+\mu} \frac{N}{N}\\
#S_1^* =& 1/R_0\\
#I_1^* =& \mu (R_0 - 1)/\beta\\
#R_1^* =& \frac{\gamma I_1^*}{(1-\Pi) + mu}\\
#R^* =& (1 - \Pi) * R_1^* / \mu\\
#Q_0=\frac{\beta_2}{\gamma+\mu} \frac{S^*}{N} + \frac{\Xi \beta_2 R_1^* I_2}{N (\gamma+\mu)}


# This creates the 'behind the scenes' code (Server)
server = function(input, output) {
  twostrain=function(t, y, parameters){
  S=ifelse(y[1]<0, 0, y[1])
  I1=ifelse(y[2]<0, 0, y[2])
  I2=ifelse(y[3]<0, 0, y[3])
  R1=ifelse(y[4]<0, 0, y[4])
  R2=ifelse(y[5]<0, 0, y[5])
  J1=ifelse(y[6]<0, 0, y[6])
  J2=ifelse(y[7]<0, 0, y[7])
  R=ifelse(y[8]<0, 0, y[8])

 with(as.list(parameters),{
  phi=(beta1*I1+beta2*I2+Theta*(beta1*J1+beta2*J2))/N
  dS = mu*N - phi *S -mu * S
  dI1= (beta1*I1+Theta*beta1*J1)*S/N - (gamma+mu)*I1
  dI2= (beta2*I2+Theta*beta2*J2)*S/N - (gamma+mu)*I2
  dR1= Pi*gamma*I1 - (beta2*I2+Theta*beta2*J2)*Xi*R1/N - mu*R1
  dR2= Pi*gamma*I2 - (beta1*I1+Theta*beta1*J1)*Xi*R2/N - mu*R2
  dJ1= (beta1*I1+Theta*beta1*J1)*Xi*R2/N - gamma*J1 - mu*J1
  dJ2= (beta2*I2+Theta*beta2*J2)*Xi*R1/N - gamma*J2 - mu*J2
  dR=(1-Pi)*gamma*(I1+I2)+gamma*(J1+J2)-mu*R
  res=c(dS, dI1, dI2, dR1,dR2,dJ1,dJ2,dR)
  return(list(res))
})
} 

  output$table1 = renderDataTable({
    R01 = (input$beta1/365)/(1/input$Ip + input$mu/365)
    S1star=1/R01
    I1star=input$mu/365 *(R01-1)/(input$beta1/365)
    R1star=1-(S1star+I1star)
    #Rstar= (1-input$pi)*R1star/input$mu/365
    R02 = (input$beta2/365)/(1/input$Ip + input$mu/365)
    Q02 = ((input$beta2/365)*S1star)/(input$mu/365 + 1/input$Ip) + (input$beta2/365)*input$xi*R1star/(input$mu/365 + 1/input$Ip) 
    S2star=1/R02
    I2star=input$mu/365 *(R02-1)/(input$beta1/365)
    R2star=1-(S2star+I2star)
    #Rstar= 1-(S1star+I1star+S1star+R1star)
    r=log(Q02)/input$Ip
    Td=log(2)/r
    data.frame(R01=round(R01,2), S1star=round(S1star,3), I1star=round(I1star,4),R02=round(R02,2), Q02=round(Q02,2),
               S2star=round(S2star,3), I2star=round(I2star,4), Td=round(Td,1))
    #  data.frame(R01=round(R01,2), S1star=round(S1star,3), I1star=round(I1star,4),R1star=round(R1star,4), R02=round(R02,2), Q02=round(Q02,2),
    #             S2star=round(S2star,3), I2star=round(I2star,4),R2star=round(R2star,4), Td=round(Td,1))
  })
  output$plot1 <- renderPlot({
    times  = seq(0, input$T[2], by=1/100)
    paras  = c(mu = input$mu, N = 1, beta1=input$beta1, beta2=input$beta2, 
               gamma = 365/input$Ip, Theta=input$theta, Xi=input$xi, Pi=input$pi)
    #paras["beta1"] = with(as.list(paras), R01*(gamma+mu))
    #paras["beta2"] = with(as.list(paras), R02*(gamma+mu))
    start = c(S = 0.06, I1 = 0.001, I2 = 0.001, R1=0, R2=0, J1=0, J2=0, R = 0.938)
    
    out = as.data.frame(ode(start, times, twostrain, paras))
    
    sel=out$time>input$T[1]&out$time<input$T[2]
    oldpar <- par(no.readonly = TRUE) # code line i
    on.exit(par(oldpar)) # code line i + 1
    par(mar = c(5,5,2,5))
    if(input$lg==FALSE){
      plot(x=out$time[sel], y=out$I2[sel], ylab="fraction", xlab="time", type="l",
           ylim=range(out[sel, c(3,4,7,8)]), xlim=c(input$T[1], input$T[2]), log=ifelse(input$lg==TRUE, "y", ""), col="red")
      lines(x=out$time, y=out$I1, col="black")
      lines(x=out$time, y=out$J1, lty=2)
      lines(x=out$time, y=out$J2, col="red", lty=2)
      
      par(new=T)
      #plot(x=out$time, y=out$S, type="l", col="green", axes=FALSE, xlab=NA, ylab=NA 
      #    ylim=range(out[sel,2]), xlim=c(input$T[1], input$T[2]), #log=ifelse(input$lg==TRUE, "y", ""))
      #axis(side = 4, col="green")
      #mtext(side = 4, line = 4, "S", col="green")
      legend("right",
             legend=c("I1", "I2", "J1", "J2"),
             lty=c(1,1,2,2),
             col=c("black", "red", "black", "red"))
    }
    
    if(input$lg==TRUE){
      plot(x=out$time[sel], y=out$I2[sel]^(1/4), ylab="fourth root fraction", xlab="time", type="l",
           ylim=range(out[sel, c(3,4,7,8)]^(1/4)), xlim=c(input$T[1], input$T[2]), col="red")
      lines(x=out$time, y=out$I1^(1/4), col="black")
      lines(x=out$time, y=out$J1^(1/4), lty=2)
      lines(x=out$time, y=out$J2^(1/4), col="red", lty=2)
      
      par(new=T)
      #plot(x=out$time, y=out$S, type="l", col="green", axes=FALSE, xlab=NA, ylab=NA 
      #    ylim=range(out[sel,2]), xlim=c(input$T[1], input$T[2]), #log=ifelse(input$lg==TRUE, "y", ""))
      #axis(side = 4, col="green")
      #mtext(side = 4, line = 4, "S", col="green")
      legend("right",
             legend=c("I1", "I2", "J1", "J2"),
             lty=c(1,1,2,2),
             col=c("black", "red", "black", "red"))
    }
  })
  
  # output$plot2b <- renderPlot({
  #   times  = seq(0, input$T[2], by=1/200)
  #   paras  = c(mu = input$mu, N = 1, beta1=input$beta1, beta2=input$beta2, 
  #              gamma = 365/input$Ip, Theta=input$theta, Xi=input$xi, Pi=input$pi)
  #   #paras["beta1"] = with(as.list(paras), R01*(gamma+mu))
  #   #paras["beta2"] = with(as.list(paras), R02*(gamma+mu))
  #   start = c(S = 0.999, I1 = 0.001, I2 = 0.00, R1=0, R2=0, J1=0, J2=0, R = 0)
  #   
  #   out1 = as.data.frame(ode(start, times, twostrain, paras))
  #   ta=out1[out1[,1]<input$T[1],]
  #   #ta=tail(out1, 1)
  #   start2 = c(S = ta[1,2], I1 = ta[1,3], I2 = 0.001, R1=ta[1,5], R2=ta[1,6], J1=ta[1,7], J2=ta[1,8], R = ta[1,9])
  #   out2 = as.data.frame(ode(start2, times, twostrain, paras))
  #   
  #   sel=out1$time>0
  #   #sel=TRUE
  #   R01 = (input$beta1/365)/(1/input$Ip + input$mu/365)
  #   R02 = (input$beta2/365)/(1/input$Ip + input$mu/365)
  #   plot(out1$S[sel], R01*out1$S[sel], ylim=c(0, R01), type="l", xlab="S", ylab="Re")
  #   lines(out2$S[sel], R01*out2$S[sel], lty=2)
  #   lines(out2$S[sel], R02*out2$S[sel], col=2)
  #   abline(h=1)
  #   
  #   legend("right",
  #          legend=c("Re1", "Re2"),
  #          col=c("black", "red"))
  # })
  # 
  output$plot2 <- renderPlot({
    times  = seq(0, input$T[2], by=1/100)
    paras  = c(mu = input$mu, N = 1, beta1=input$beta1, beta2=input$beta2, 
               gamma = 365/input$Ip, Theta=input$theta, Xi=input$xi, Pi=input$pi)
    #paras["beta1"] = with(as.list(paras), R01*(gamma+mu))
    #paras["beta2"] = with(as.list(paras), R02*(gamma+mu))
    start = c(S = 0.999, I1 = 0.001, I2 = 0.00, R1=0, R2=0, J1=0, J2=0, R = 0)
    
    out1 = as.data.frame(ode(start, times, twostrain, paras))
    ta=out1[out1[,1]<input$T[1],]
    #ta=tail(out1, 1)
    start2 = c(S = ta[1,2], I1 = ta[1,3], I2 = 0.001, R1=ta[1,5], R2=ta[1,6], J1=ta[1,7], J2=ta[1,8], R = ta[1,9])
    out2 = as.data.frame(ode(start2, times, twostrain, paras))
    
    sel=out1$time>0
    #sel=TRUE
    R01 = (input$beta1/365)/(1/input$Ip + input$mu/365)
    R02 = (input$beta2/365)/(1/input$Ip + input$mu/365)
    plot(out1$S[sel], R01*out1$S[sel], ylim=c(0, R01), type="l", xlab="S", ylab="Re")
    lines(out2$S[sel], R01*out2$S[sel], lty=2)
    lines(out2$S[sel], R02*out2$S[sel], col=2)
    abline(h=1)
    points(1/R01, 1, pch="X")
    points(1/R02, 1, pch="O")
    
    legend("right",
           legend=c("Re1", "Re2"),
           lty=c(1,1),
           col=c("black", "red"))
  })
  
  output$plot3 <- renderPlot({
    times  = seq(0, input$T[2], by=1/200)
    paras  = c(mu = input$mu, N = 1, beta1=input$beta1, beta2=input$beta2, 
               gamma = 365/input$Ip, Theta=input$theta, Xi=input$xi, Pi=input$pi)
    #paras["beta1"] = with(as.list(paras), R01*(gamma+mu))
    #paras["beta2"] = with(as.list(paras), R02*(gamma+mu))
    start = c(S = 0.999, I1 = 0.001, I2 = 0.00, R1=0, R2=0, J1=0, J2=0, R = 0)
    
    out1 = as.data.frame(ode(start, times, twostrain, paras))
    ta=out1[out1[,1]<input$T[1],]
    #ta=tail(out1, 1)
    start2 = c(S = ta[1,2], I1 = ta[1,3], I2 = 0.001, R1=ta[1,5], R2=ta[1,6], J1=ta[1,7], J2=ta[1,8], R = ta[1,9])
    out2 = as.data.frame(ode(start2, times, twostrain, paras))
    
    sel=out1$time>0
    #sel=TRUE
    R01 = (input$beta1/365)/(1/input$Ip + input$mu/365)
    R02 = (input$beta2/365)/(1/input$Ip + input$mu/365)
    plot(out2$S[sel], out1$I1[sel], log="", xlim=c(0, 1), type="l", xlab="S", ylab="I")
    lines(out2$S[sel], out2$I1[sel], lty=2)
    lines(out2$S[sel], out2$I2[sel], col=2)
    # abline(h=1)
    
    legend("right",
           legend=c("I1only", "I2", "I1tostrain"),
           lty=c(1,1,2),
           col=c("black", "red", "black"))
  })
}
)