[R-sig-dyn-mod] iterative procedure vs deSolve
Paula Sanchez Marin
paulasanchez at uvigo.es
Thu Mar 17 14:56:40 CET 2016
Thank you Daniel,
The analytical solution of the ODE is a system of two non-linear
equations with two variables.
You say I can use those expressions and evaluate them at my time points
with the parameter values and compare the result with the output from ode()
That's what I was trying to do in excel. The way I know of solving this
is by using an iterative procedure. If there is another way I am missing
I will be glad to know!
Since Crest is needed to calculate Cgills and Cgills is needed to
calculate Crest, there is a "circular reference" there, that excel
solves by performing a repeated recalculation of the worksheet cells
containing the formulas until a specific numeric condition is met (until
the change is minimal).
Cheers,
Paula
> ------------------------------
>
> 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.
>>
>> _______________________________________________
>> R-sig-dynamic-models mailing list
>> R-sig-dynamic-models at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models
>
>
> ------------------------------
>
> Subject: Digest Footer
>
> _______________________________________________
> R-sig-dynamic-models mailing list
> R-sig-dynamic-models at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models
>
> ------------------------------
>
> End of R-sig-dynamic-models Digest, Vol 101, Issue 1
> ****************************************************
--
Paula Sánchez Marín
Dep. Ecoloxía e Bioloxía Animal
Laboratorio Torre CACTI 99
Edificio Citexvi
Universidade de Vigo
Tel: 986818780
mailto: paulasanchez at uvigo.es
More information about the R-sig-dynamic-models
mailing list