Bioequivalence and Bioavailability Forum • 1 compartment model with lag time

1 compartment model with lag time [PK / PD]

posted by yjlee168 – Kaohsiung, Taiwan, 2013-07-17 22:57 (4711 d 03:17 ago) – Posting: # 11005
Views: 25,263
(edited on 2013-07-18 13:00)

Hi John
Let me try to solve this puzzle. I used Helmut's equation since his is correct one. What I am using is JAGS with R. Two R packages will be required (rjags and R2jags). Here is its brief tutorial. Many can be found with Google. Since you don't give any value for this question, I use following parameters to do simulation dataset first.

Parameter values I used:

D = 100

V = 10

K01 = 1.67  <-- we want to calculate this param.

K10 = 0.213

t_lag = 0.5

So I got C(t) as follows:


t (hr)  C(t)

0       0

1       5.331

2       7.391

3       6.553

4       5.405

6       3.551 <- we will use this C(t).

8       2.320

12      0.990

16      0.422

18      0.276

24      0.0768

Here is the JAGS model. You can copy & paste the model and save it as "jag009.txt" in your R's working directory (or folder).


model

    {

     C~dnorm(mu, 1.0E+6)     # C: drug concentration

     mu<-D*K01/(V*(K01-K10))*(exp(-K10*(t-t_lag))-exp(-K01*(t-t_lag)))

     K01~dnorm(theta1, 10)   # define K10 = theta1 here

     theta1<-0.8             # prior of K01; as close to the real number

                             # as possible; or just guess one if don't know it.

    }

Then run the following R script.


library(rjag)

library(R2jags)

dataList= list(

C=3.551,

D=100,

V=10,

K10=0.213,

t=6,

t_lag=0.5

)

params = c("K01")              # The parameter(s) to be monitored.

initsList = list(K01=0.8)      # initialize prior here;I tried not to give too close value for K01

adaptSteps = 500               # Number of steps to "tune" the samplers.

burnInSteps = 6000             # Number of steps to "burn-in" the samplers.

nChains = 3                    # Number of chains to run.

numSavedSteps=100000           # Total number of steps in chains to save.

thinSteps=1                    # Number of steps to "thin" (1=keep every step).

nIter = ceiling(( numSavedSteps * thinSteps ) / nChains )   # Steps per chain.

# Create, initialize, and adapt the model:

jagsModel = jags.model("jag009.txt", data=dataList ,    

n.chains=nChains , n.adapt=adaptSteps, inits=initsList) 

# Burn-in:

cat("Burning in the MCMC chain...\n")

update(jagsModel, n.iter=burnInSteps)

# The saved MCMC chain:

cat("Sampling final MCMC chain...\n")

codaSamples <- coda.samples(jagsModel, params, n.iter=nIter)

show(summary(codaSamples))

dev.new()

plot(codaSamples)

dev.new()

gelman.plot(codaSamples)

cat("\n")

K01.mat <- as.matrix(codaSamples[[1]])

posterior.estimates <- rbind(K01.mat)

vars <- t(as.matrix(posterior.estimates))

params.mean <- apply(vars, 1, mean)

K01.mean <- as.matrix(params.mean[[1]])

C<-3.551; D<-100; V<-10; K10<-0.213; t<-6; t_lag<-0.5  ### all known parameter values here

C_calc<-D*K01.mean/(V*(K01.mean-K10))*(exp(-K10*(t-t_lag))-exp(-K01.mean*(t-t_lag)))

out<-data.frame(Parameters=c("sim. C(t)","calc. C(t)"),Values=c(C, C_calc))

cat("\n");show(out);cat("\n\n")

Finally, I got


...

Iterations = 6501:39834

Thinning interval = 1 

Number of chains = 3 

Sample size per chain = 33334 



1. Empirical mean and standard deviation for each variable,

   plus standard error of the mean:



          Mean             SD       Naive SE Time-series SE

     1.670e+00     3.281e-03      1.038e-05      1.319e-05



2. Quantiles for each variable:



 2.5%   25%   50%   75% 97.5%

1.663 1.667 1.670 1.672 1.676



  Parameters   Values

1  sim. C(t) 3.551000

2 calc. C(t) 3.551006

It converges successfully even with init. K01 = 0.8 from two diagnostic plots (not attached with this post). Try it yourself. Good luck.

—
All the best,
-- Yung-jin Lee
bear v2.9.6:- created by Hsin-ya Lee & Yung-jin Lee
Kaohsiung, Taiwan https://www.pkpd168.com/bear
Download link (updated) -> here

Complete thread:

1 compartment model with lag time jag009 2013-07-17 15:58 [PK / PD]
- solve for k01? Helmut 2013-07-17 16:08
  - solve for k01? jag009 2013-07-17 16:19
    - solve for k01? ElMaestro 2013-07-17 17:09
    - no closed form Helmut 2013-07-17 23:01
- 1 compartment model with lag timeyjlee168 2013-07-17 20:57
- Package-free solution ElMaestro 2013-07-17 21:53
  - Package-free solution yjlee168 2013-07-17 22:11
    - Package-free solution ElMaestro 2013-07-17 22:27
      - Package-free solution yjlee168 2013-07-17 22:42
  - Package-free solution ElMaestro 2013-07-17 23:21
    - A new R-King was born d_labes 2013-07-18 08:45
      - Brent ElMaestro 2013-07-18 09:07
      - A new R-King was born yjlee168 2013-07-18 10:21
        
        OT: TTT in bear Helmut 2013-07-18 20:28
        
        OT: TTT in bear yjlee168 2013-07-18 22:17
        
        OT: TTT subtleties d_labes 2013-07-19 09:47
        
        OT: TTT subtleties yjlee168 2013-07-20 19:44
        
        OT: beyond TTT? Helmut 2013-07-21 00:08
        
        OT: beyond TTT? yjlee168 2013-07-22 13:45
        
        OT: keep TTT! Helmut 2013-07-22 14:21
        
        OT: keep TTT! yjlee168 2013-07-22 23:42
        
        OT: EOD (here) Helmut 2013-07-23 01:06
    - Package-free solution yjlee168 2013-07-18 10:27
      - Thank you guys! jag009 2013-07-18 16:37
        
        Thank you guys! ElMaestro 2013-07-18 18:57
        
        Thank you guys! jag009 2013-07-19 20:40