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

[Daniel Kaschek]

Dear Paula,

I did not really understand what you mean by "iterative procedure" in 
Excel. You write that it is based on the analytical solution of the ODE:

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))

Have you checked the solution? You could use your expressions above and 
evaluate them at your time points with the parameter values and compare 
the result with the output from ode(). If they do not coincide, either 
your analytical solution is wrong or you misspecified the ODE.

Cheers,
Daniel


On Mi, Mär 16, 2016 at 7:59 , Paula Sanchez Marin 
<paulasanchez at uvigo.es> wrote:
> 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.
> 
> _______________________________________________
> R-sig-dynamic-models mailing list
> R-sig-dynamic-models at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models