model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    B <- N - S - D - L
    dtS <- iS*B - eS*S
    dtD <- fi*iL*B - fe*eL*D 
    dtL <- (nLN-1)*iL*B - eL*L
    return(list(c(dtS, dtD, dtL)))  
  }) 
}  

# B: blood, S: spleen, D: draining LN, L: LN
# Parameters based upon Textor et al PLOS Comp Biol 2014

nLN <- 39
p <- c(N=100,iS=1,iL=1.5/nLN,eS=1/6,eL=1/13.5,fi=1,fe=1)
s <- c(S=0,D=0,L=0)
tweak <- "nsol$B=parms[1]-nsol$S-nsol$D-nsol$L"
s <- run(tweak=tweak)
s <- newton(s)
p["fe"] <- 0
run(24*1,tweak=tweak)
p["fi"] <- 9
run(24*1,tweak=tweak)