------------------------------
Message: 2
Date: Thu, 17 Mar 2016 10:24:54 +0100
From: Daniel Kaschek <daniel.kaschek at physik.uni-freiburg.de>
To: Special Interest Group for Dynamic Simulation Models in R
<r-sig-dynamic-models at r-project.org>
Subject: Re: [R-sig-dyn-mod] iterative procedure vs deSolve
Message-ID: <1458206694.5529.0 at physik.uni-freiburg.de>
Content-Type: text/plain; charset=utf-8; format=flowed
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.