[R] Fitting compartmental model with nls and lsoda?
DivineSAAM@aol.com
DivineSAAM at aol.com
Fri Jan 23 00:56:02 CET 2004
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
More information about the R-help
mailing list