model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    dN <- (pN-d)*N
    return(list(dN))  
  }) 
}  

sm <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    tlag <- t - Delta
    if (tlag < 0) lags <- 0
    else lags <- lagvalue(tlag,1)    # return lag of A
    dA <- -(p+d)*A + 2*p*lags[1]*exp(-Delta*d)
    dB <- p*A - d*B -  p*lags[1]*exp(-Delta*d)
    return(list(c(dA, dB)))  
  }) 
}  

erl <- function(t, state, parms){
  with(as.list(c(state,parms)),{
    gamma <- ndiv/Delta
    B <- state[1:ndiv]
    dB <- rep(0,ndiv)
    dB[1] <- p*A - (d+gamma)*B[1]
    dB[2:ndiv] <- gamma*B[1:(ndiv-1)] - (d+gamma)*B[2:ndiv]
    dA <- -(p+d)*A + 2*gamma*B[ndiv]
    return(list(c(dB,dA)))    
  }) 
} 

# Ganusov et al. J. Immunol. Methods 2005:
# Solve net growth rate (r) of SM model from 0 = 2pe^[-(d+r)Delta] - (p+d+r)

rsm <- function(r,parms){
  with(as.list(parms),return(2*p*exp(-(d+r)*Delta) - (p+d+r)))
}

p <- c(p=1,d=0.01,Delta=0.5,pN=0)
uniroot.all(rsm,c(0,10),parms=p)    # growth rate
p["pN"] <- 1/(p["Delta"]+1/p["p"])  # equivalent division rate in ODE

s <- c(N=1)      # state of ODE
run(20,tstep=1/24,log="y")

s <- c(A=1,B=0)  # state of Smith-Martin model
run(20,tstep=1/24,odes=sm,delay=TRUE,tweak="nsol$N=nsol$A+nsol$B",add=TRUE,show="N")

ndiv <- 10
B <- rep(0,ndiv)
names(B) <- paste("B",seq(1,ndiv),sep="")
s <- c(B,A=1)    # state of Erlang model
run(20,tstep=1/24,odes=erl,tweak="nsol$N=nsol$A+rowSums(nsol[,2:(ndiv+1)])",add=TRUE,show="N")