R help - Jan 2004 - Fitting compartmental model with nls and lsoda?

Dear Colleagues,

Our group is also working on implementing the use of R for pharmacokinetic
compartmental analysis. Perhaps I have missed something, but
> fit <- nls(noisy ~ lsoda(xstart, time, one.compartment.model, c(K1=0.5,
k2=0.5)),
+            data=C1.lsoda,
+            start=list(K1=0.3, k2=0.7),
+            trace=T
+            )
Error in eval(as.name(varName), data) : Object "C1.lsoda" not found

What part of the e-mail did I miss? I would like to get this problem up an
running.

Now, I am including Richar Upton's 2 cm model implementation and Christoffer
Tornoe's nls solution (I recommend Christoffer's nlmeODE package for
these problems also if multi-response data is available) The code follows:
--------------------------------------------------------------
# Simulation of a 2 compartment pharmacokinetic model using "R"
# Richard N. Upton, 11/3/02, richard.upton at adelaide.edu.au

# The "R" home page is at http://www.R-project.org/

# I make no representations about being an "R" guru.  My contribution
here
# is hopefully to provide a starting point in "R" for people who
# have a pharmacokinetic modelling background.

# This text can  be "cut and pasted" into R, or read in as a
"source" file

# There are two differential equations in the system:

# V*dC/dt = Doserate - C*Cl + k21*A2 - k12*V*C
# dA2/dt = k12*V*C - k21*A2

# C is a dependent variable (Concentration in the central compartment)
# A is a dependent variable (Amount in the second compartment)
# t is the independent variable (time)

# V is the volume of the central compartment
# Cl is the clearance from the central compartment
# k12 is the rate constant between the central and second compartment
# k21 is the rate constant between the second and central compartment

# Dose is the total amount of drug given
# tau is the time over which this amount is given
# The doserate (amount/time) is therefore Dose/tau
# A bolus dose should be thought of as a short infusion

# The lsoda function is very fussy about names for variables
# C[1] = C, meaning the first dependent variable ; Cd1 is its derivative wrt
time
# C[2] = A2, meaning the second dependent variable ; Cd2 is its derivative wrt
time
# You can change C to another name, but must keep these conventions
# The output from Cprime (its last line) must be a list of the derivative of C
wrt time

# You must install the "odesolve" package in R.  See the website for
details.

# This example gave results similar (within 6 sig. fig.) to the same problem
# solved in an independent differential equation solving package.


#Load the odesolve package
require(odesolve)

#Specify parameters
times <- c(0:180)
tau  <- 4
Dose <- 30
V    <- 23.1
Cl   <- 1
k12  <- 0.197118
k21  <- 0.022665

#A quick check - compare these steady-state values with values after a long
infusion
Css <- (Dose/tau)/Cl
A2ss <- V*Css*(k12/k21)

#lsoda requires the parameters as an object (p) with names
p <- c(V=V, Cl=Cl, k12=k12, k21=k21)


#Differential equations are declared in a function
Cprime <- function(t, C, p)
 
          {
          if (t < tau) Doserate <- (Dose/tau) else Doserate <- 0
          Cd1 <- (Doserate - C[1]*p["Cl"] + p["k21"]*C[2]
- p["V"]*p["k12"]*C[1])/p["V"]
          Cd2 <- (p["V"]*p["k12"]*C[1] -
p["k21"]*C[2])
          list(c(Cd1, Cd2))
          }
      
 
#Solve the system of differential equations, including initial values
result <- lsoda( c(0,0), times, Cprime, p)

#Reformat the result
result <- data.frame(result)
names(result) <- c("time","Conc", "Amount2")

#Have a look at the result
print(result)
plot(result$Conc ~ result$time, type="b", main="Central
compartment", xlab="time", ylab="Conc")
plot(result$Amount2 ~ result$time, type="b", main="Second
compartment", xlab = "time", ylab = "Amount")
--------------------------------------------------------------

Our group is also working on implementing a ODE solvers suite for R for small to
medium size problems.

Thanks! for bringing this type of discussion to the R-news.

olinares at med.umich.edu

Oscar A. Linares, MD, PhD              ///////
Michigan Diabetes Institute S c I S O F T ///=20=03
Molecular Medicine Unit       ______////////=20=04
SciSoft Group                  \_\_\_\/////
Ann Arbor, MI                   \_\_\_\///
Tel. (734) 637-7997              \_\_\_\/
Fax. (734) 453-7019

> Now, I am including Richar Upton's 2 cm model implementation and
> Christoffer Tornoe's nls solution (I recommend Christoffer's
> nlmeODE package for these problems also if multi-response data is
> available) The code follows:
Multi-response data?, Is there a way of dealing with those multiresponse
nonlinear regression problems in R (a la Chapter 4 of Bates and Watts, if
possible)?, the only alternative I heard of was in a mail from Douglas
Bates, using optim() and the definition of the Box-Drapper criterion,

 prod(svd(y-f(p1,p2,p3), nu=0, nv=0)$d)^2

if you wish to minimize a matrix of y's in respect to p1,p2 and p3.


best regards,

Jesus

--------------------------------------------------------------
Jes?s Mar?a Fr?as Celayeta
School of Food Sci. and Env. Health.
Faculty of Tourism and Food
Dublin Institute of Technology
Cathal Brugha St., Dublin 1. Ireland
Phone: +353 1 4024459 Fax: +353 1 4024495
http://www.dit.ie/DIT/tourismfood/science/staff/frias.html
--------------------------------------------------------------


-- 
This message has been scanned for content and 
viruses by the DIT Information Services MailScanner 
Service, and is believed to be clean.
http://www.dit.ie

Dear Jesus,

Yes there is a way and it is via Christoffer Torn\{o}e's package nlmeODE. I
checked Chapter 4 of Bates and Watts. As usual, a little work involved, but
doable and quite powerful.

Thanks,
olinares

Dear Jesus

It is possible with nlmeODE. Check out the example ?PKPDpool where PK and PD
are fitted simultaneously.

Best regards

Christoffer
> Now, I am including Richar Upton's 2 cm model implementation and 
> Christoffer Tornoe's nls solution (I recommend Christoffer's
nlmeODE
> package for these problems also if multi-response data is
> available) The code follows:
>Multi-response data?, Is there a way of dealing with those multiresponse
nonlinear regression problems in R (a la Chapter 4 of Bates and Watts,
if
>>possible)?, the only alternative I heard of was in a mail from Douglas
Bates, using optim() and the definition of the Box-Drapper criterion,
> prod(svd(y-f(p1,p2,p3), nu=0, nv=0)$d)^2
>if you wish to minimize a matrix of y's in respect to p1,p2 and p3.

best regards,

Jesus


**********************************************************************
Proprietary or confidential information belonging to Ferring Holding SA or to
one of its affiliated companies may be contained in the message. If you are not
the addressee indicated in this message (or responsible for the delivery of the
message to such person), please do not copy or deliver this message to anyone.
In such case, please destroy this message and notify the sender by reply e-mail.
Please advise the sender immediately if you or your employer do not consent to
e-mail for messages of this kind. Opinions, conclusions and other information in
this message represent the opinion of the sender and do not necessarily
represent or reflect the views and opinions of Ferring.
**********************************************************************


	[[alternative HTML version deleted]]

>>>>> "DivineSAAM" == DivineSAAM  <DivineSAAM at
aol.com> writes:
    > Dear Colleagues,
    > Our group is also working on implementing the use of R for
    > pharmacokinetic compartmental analysis. Perhaps I have
    > missed something, but 

    > > fit <- nls(noisy ~ lsoda(xstart, time,
    > one.compartment.model, c(K1=0.5, k2=0.5)), 
    > +            data=C1.lsoda,
    > +            start=list(K1=0.3, k2=0.7),
    > +            trace=T
    > +            )
    > Error in eval(as.name(varName), data) : Object "C1.lsoda" not
found

    > What part of the e-mail did I miss? I would like to get
    > this problem up an running. 

Oscar,

There were several problems with the code in my previous posting.
I'll append an example that does work.  While this example often
works, there are cases when nls fails by reducing the step factor
below minFactor.  It is even less stable for fitting less trivial
examples with real data sets.  I'll be very interested to keep in
touch as we make progress on this problem.

Regards, Mike


Here's my (more) correct example:

require(odesolve)

## Simple one compartment blood flow model:
one.compartment.model <- function(t, x, parms) {
  C1 <- x[1] # compartment
  with(as.list(parms),{
    input <- approx(signal$time, signal$input, t)$y
    dC1 <- K1 * input - k2 * C1
    list(c(dC1))
  })
}

## vector of timesteps
time <- seq(0, 100, 1)

## external signal with rectangle impulse
signal <- as.data.frame(list(time=time,
                             input=rep(0,length(time))))
signal$input[signal$time >= 10 & signal$time <=40] <- 0.2

## Parameters for steady state conditions
parms <- c(K1=0.5, k2=0.5)

## Start values for steady state
xstart <- c(C1=0)

## calculate C1 with lsoda:
C1.lsoda <- as.data.frame(lsoda(xstart, time, one.compartment.model, parms))

## Add some noise to the output curve
C1.lsoda$noisy <- C1.lsoda$C1 + rnorm(nrow(C1.lsoda), sd=0.15*C1.lsoda$C1)

## See if I can run a fit to find the parameters that I started with...
require(nls)
fit <- nls(noisy ~ lsoda(xstart, time, one.compartment.model, c(K1=K1,
k2=k2))[,2],
           data=C1.lsoda,
           start=list(K1=0.3, k2=0.7),
           trace=T,
           control=list(tol=1e-2,
             minFactor=1/1024/1024)
           )

fit.rk4 <- nls(noisy ~ rk4(xstart, time, one.compartment.model, c(K1=K1,
k2=k2))[,2],
               data=C1.lsoda,
               start=list(K1=0.3, k2=0.7),
               trace=T,
               control=list(tol=1e-2,
                 minFactor=1/1024/1024)
               )

## Plot what I've got so far:
par(mfrow=c(2,2))
plot(noisy ~ time, data=C1.lsoda, main='Input, C1, C1+noise',
col='forestgreen')
points(input ~ time, data=signal, type='b')
points(C1 ~ time, data=C1.lsoda, type='b', pch=16)

t <- seq(0,100,0.1)
plot(noisy ~ time, data=C1.lsoda, main='Input, C1+noise, lsoda',
     col='forestgreen')
lines(t, predict(fit, list(time=t)), col='red',type='l')

plot(noisy ~ time, data=C1.lsoda, main='Input, C1+noise, rk4',
     col='forestgreen')
lines(t, predict(fit.rk4, list(time=t)), col='blue', type='l')

plot(t, predict(fit.rk4, list(time=t))-predict(fit, list(time=t)),
     type='l')
abline(h=0)


print('Coefficients for lsoda solution:')
print(coef(fit))
print(vcov(fit))

print('Coefficients for rk4 solution:')
print(coef(fit.rk4))
print(vcov(fit.rk4))

-- 
Michael A. Miller                               mmiller3 at iupui.edu
  Imaging Sciences, Department of Radiology, IU School of Medicine