model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    dtT <- s - dT*T - beta*(1-epsB)*T*V
    dtI1 <- beta*(1-epsB)*T*V - (d1 + gamma*(1-epsG) + k1)*I1
    dtI2 <- gamma*(1-epsG)*I1 - (d2 + k2)*I2
    dtV <- c*I2 - c*V
    return(list(c(dtT, dtI1, dtI2, dtV)))  
  }) 
}  

# Parameters based upon Gadhamsetty et al. J. Virol. 2016 (late killing regime)
# But not making the QSSA V'=0 and setting p=c=25 to scale V to I2.
# See Cardoza et al. PLOS Pathogens 2017 for experimental data

# Parameters for late killing:
p <- c(s=0.1,dT=0.1,d1=0.1,d2=2,k1=0,k2=0,beta=9.1,gamma=1,c=24,epsB=0,epsG=0)
s <- c(T=1,I1=0,I2=0,V=0)
with(as.list(p),(d1+gamma+1.5)*(d2+1.5)/gamma)
newton(s)
continue(s,x="beta",y="V",ymin=-0.01)
p["k2"] <- 3
newton(s)

s <- c(T=1,I1=0,I2=0,V=1)
chronic <- newton(run(100))

# Start conventional therapy in chronic state (90% effective):
p["epsB"] <- 0.9; f1 <- run(14,state=chronic,tstep=1/24,ymin=1e-4,log="y",show="V")
log(f1["V"]/chronic["V"])/14         # Check if downslope is close to 1/day
newton(chronic)
# Add an intergrase inhibitor to the same therapy:
p["epsG"] <- 0.9; f2 <- run(14,state=chronic,tstep=1/24,show="V",add=TRUE)
newton(chronic)

# Parameters for early killing:
p <- c(s=0.1,dT=0.1,d1=1,d2=1,k1=0,k2=0,beta=8.75,gamma=1,c=24,epsB=0,epsG=0)
with(as.list(p),(d1+gamma+1.5)*(d2+1.5)/gamma)
p["k1"] <- 3
chronic <- newton(run(100))

# Start conventional therapy in chronic state (90% effective):
p["epsB"] <- 0.9; f1 <- run(14,state=chronic,tstep=1/24,ymin=1e-4,log="y",show="V")
log(f1["V"]/chronic["V"])/14         # Check if downslope is close to 1/day
newton(chronic)
# Add an intergrase inhibitor to the same therapy:
p["epsG"] <- 0.9; f2 <- run(14,state=chronic,tstep=1/24,show="V",add=TRUE)
newton(chronic)