[R-sig-dyn-mod] iterative procedure vs deSolve

[Paula Sanchez Marin]

Hello,

While trying to fit some experimental data to a pharmacokinetic model I 
have observed that the results of the model predicted by deSolve are not 
the same as those resulting from plotting of the solved equations using 
an iterative procedure.
I would like to know why do the two methods give different results and 
which one is the correct one.

Here it goes my example:
A contaminant (Cu) enter an organism by compartment 1 "gills" and from 
there it can be transferred to compartment 2 "rest".
The contaminant can only be eliminated from compartment 1 to the 
environment.
The parameters (transfer coefficients) are: k01, k10, k12 and k21.
I introduced the concentration or dilution factors when Cu is transfered 
from 1 to 2 and viceversa by using the ratios between the weight of the 
gills (wgills) and the weight of the rest (wrest).

The differential equations describing changes in concentration in both 
compartments with time are:
dCgills/dt = k01*Cu+k21*wrest/wgills*Crest-(k10+k12)*Cgills
dCrest/dt=k12*wgills/wrest*Cgills-k21*Crest

This is the script I used for solving the differential equations in R:

#Constants
Cu=0.03
wrest=0.18
wgills=0.015
#Model
mod<-function(t,y,parms){
   with(as.list(c(parms,y)), {
     dCgills<-k01*Cu+k21*wrest/wgills*Crest-(k10+k12)*Cgills
     dCrest<-k12*wgills/wrest*Cgills-k21*Crest
     list(c(dCgills,dCrest))
   })
}
#Parameters
parms=c(k10=0.211,k12=0.234,k01=1196,k21=0.073)
#Initial conditions
X0<-c(Cgills=5.771,Crest=5.905)
#Output
out<-as.data.frame(lsoda(y=X0,times=0:100,func=mod,parms))
write.csv2(out,"out.csv")

I used the same Constants (Cu, wrest and wgills), and the same 
Parameters and Initial conditions and calculated in excel (allowing 
iterations) the results using the equations below (the solution to the 
differential equations):
Cgills=Cgills_ini*exp(-(k01+k12)*t)+(Cu*k01+Crest*k21*wrest/wgills)/(k10+k12)*(1-exp(-(k10+k12)*t))
Crest=Crest_ini*exp(-k21*t)+Cgills*k12*wgills/wrest/k21*(1-exp(-k21*t))


Here are the results (for the first 10 hours) from R:

     time    Cgills     Crest
1      0   5.77100  5.905000
2      1  36.80185  5.914667
3      2  56.88706  6.394786
4      3  70.21069  7.149393
5      4  79.34640  8.058456
6      5  85.87763  9.048616
7      6  90.77894 10.075197
8      7  94.65042 11.111150
9      8  97.86195 12.140248
10     9 100.64178 13.152906
11    10 103.13106 14.143612

And from the iterative procedure in excel:
t 	Cgills 	Crest
0 	5,771 	5,905
1 	37,0318579 	6,18568653
2 	58,2609877 	7,21694678
3 	73,4040942 	8,60007128
4 	84,5860141 	10,1314063
5 	93,0882873 	11,7033524
6 	99,7407581 	13,260262
7 	105,102356 	14,7754342
8 	109,556482 	16,2380253
9 	113,367627 	17,6454062
10 	116,718031 	18,998768


And here the plot:


Why do the curves differ?

Thank you very much for your guindance.