Population values not updating in deSolve in R

Question

I am working on a model with different compartments for different age groups and smokers/non-smokers. I thought all was going well, but when I started to dig into the output, I see that the populations are not updating at each time step as I anticipated. I used the debugger to step through the code, and indeed, it keeps using the initial stock size for every time step, so I end up with a constant change in population sizes, which is incorrect.

There are 200 stocks in the model, so I don't want to paste the whole thing in, but here are the pieces I think may be most relevant.

smokdiff <- function(times,init,parameters) {
  with(as.list(c(init,parameters)), {

    dns00 <- births - ns00*mort[1,1] - ns00/5 
    dns05 <- ns00/5 - ns05*mort[1,2] - ns05/5

...

    list(c(dns00,dns05,dns10,...
 })
}

times <- seq(0,tf,dt)
init <- c(ns00,ns05,ns10,...
parameters <- c(births,mort,init15,utinit15,ltinit15)
out <- ode(func=smokdiff, times=times, y=init, parms=parameters)

Here is a sample of the output for the ns00 population:

     time       1
[1,]    0 9244769
[2,]    1 9217034
[3,]    2 9189300
[4,]    3 9161566
[5,]    4 9133832
[6,]    5 9106097 

As you can see, it decreases by ~27,734 each step, which is only correct for the first step (e.g. it's supposed to decrease by 22,104 in the second step).

I have compared to multiple models I've found online, and I can't spot what it is I'm doing wrong. I've also tried other "methods" in the ode call, but they all display the constant change. I'm beginning to worry that I have some basic misunderstanding about how this is supposed to work. Can anyone see what I'm doing wrong? Thank you!

SOLUTION: The initial states have to be defined in a 'state=value' manner. Simply putting a variable with the value in it is insufficient. Here is a small working example:

library(deSolve)
smokdiff <- function(times,init,parameters) {
  with(as.list(c(init,parameters)), {

    dns00 <- births - ns00*mort[1] - ns00/5 
    dns05 <- ns00/5 - ns05*mort[2] - ns05/5
    dns10 <- ns05/5 - ns10*mort[3] - ns10/5

    list(c(dns00,dns05,dns10))
 })
}

mort <- c(0.003, 0.001, 0.001)
ns00 <- 9000000; ns05 <- 8500000; ns10 <- 8500000
births <- ns00/5
times <- seq(0,50,0.25)
init <- c(ns00=ns00,ns05=ns05,ns10=ns10)
parameters <- c(births,mort)
out <- ode(func=smokdiff, times=times, y=init, parms=parameters)

Answer 1

A note for any passing readers - this comes down to how ode() and with() interact.

with() goes looking for variable names, and if they are not present, fails silently.

That's why the solution above works, it gives with() those variable names, but this comes down to poor design in the with() function really.