[R] Fitting compartmental model with nls and lsoda?

[DivineSAAM@aol.com]

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