IV infusion with ODE models

[Sibojang9]

Hi Andrew,

It seems that how to implement IV infusion with ODE models is not documented in PopEd or Desolve packages. Could you please shed some light on that?

Thank you,

Sibo

[JeanUTHSC]

I agree, iv infusion is a highly desirable function. It's going to be really helpful to have it in PopED :)

[andrewhooker]

Hi

This is certainly possible. I will try to send an example in the next day or two using deSolve. You can also use mrgsolve, pkpdsim or rxode to define your ode model including iv infusion and use that to drive poped if you like.

Best regards, Andrew

Andrew Hooker, Ph.D. Associate Professor of Pharmacometrics Dept. of Pharmaceutical Biosciences Uppsala University Box 591, 751 24, Uppsala, Sweden Phone: +46 18 471 4355 Mobile: +46 768 000 725 www.farmbio.uu.se/research/researchgroups/pharmacometrics/

On 8 Jul 2017, 18:04 +0200, JeanUTHSC notifications@github.com, wrote:

I agree, iv infusion is a highly desirable function. It's going to be really helpful to have it in PopED :) — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub, or mute the thread.

[Sibojang9]

Hi Andrew,

Thank you so much for providing more Desolve examples. Pkpdsim and Rxode are definitely great alternatives since they are specially designed for pharmacometrics and already widely used in this community. The thing is that I am not very clear how PopDE communicates with those tools, i.e. how PopDE feeds the ODE and receives the feedback. At this time point, I think that deSolve would fit my project well. Use of PopDE in conjunction with Pkpdsim or Rxode would be a fantastic feature of PopDE. I would be excited to see more tutorials about this later on. I really appreciate that you spearhead this area.

Sibo

[andrewhooker]

Hi

Below is an example with IV infusion.

library(PopED)

library(deSolve)

PK_1comp_ode_inf <- function(Time, State, Pars){
  with(as.list(c(State, Pars)), {
    if( Time %% TAU < TINF ){ 
      In  <-  DOSE/TINF
    } else {
      In  <-  0
    }
    dA1 <- In - CL/V*A1
    return(list(c(dA1)))
  })
}

PK_1comp_ode_inf_ff <- function(model_switch, xt, parameters, poped.db){
  with(as.list(parameters),{
    A_ini <- c(A1=0)

    times_xt <- drop(xt)
    times <- c(0,times_xt) ## add extra time for start of integration
    times <- sort(times) 
    times <- unique(times) # remove duplicates

    out <- ode(A_ini, times, PK_1comp_ode_inf, parameters)#, events = list(data = eventdat))#atol=1e-13,rtol=1e-13)
    y = out[, "A1"]/(V)
    y=y[match(times_xt,out[,"time"])]
    y=cbind(y)
    return(list(y=y,poped.db=poped.db))
  })
}

## -- parameter definition function 
## -- names match parameters in function ff
sfg <- function(x,a,bpop,b,bocc){
  parameters=c( CL=bpop[1]*exp(b[1]),
                V=bpop[2]*exp(b[2]),
                DOSE=a[1],
                TINF=a[2],
                TAU=a[3])
  return( parameters ) 
}

poped.db <- create.poped.database(ff_file = "PK_1comp_ode_inf_ff",
                                  fError_file="feps.add.prop",
                                  fg_file="sfg",
                                  groupsize=20,
                                  m=1,      #number of groups
                                  sigma=c(0.04,0.01),
                                  notfixed_sigma = c(1,0),
                                  bpop=c(CL=20,V=72.8), 
                                  d=c(CL=0.09,V=0.09), 
                                  xt=c(1,2,8,12,24),
                                  minxt=0,
                                  maxxt=100,
                                  a=c(DOSE=100,TINF=10,TAU=24))

plot_model_prediction(poped.db)

# slow evaluation time
tic(); evaluate_design(poped.db); toc()

##################################
# Compiled ODE
##################################
library(Rcpp)

cppFunction('List one_comp_ode(double Time, NumericVector A, NumericVector Pars) {
            int n = A.size();
            NumericVector dA(n);

            double CL = Pars[0];
            double V = Pars[1];
            double DOSE = Pars[2];
            double TINF = Pars[3];
            double TAU = Pars[4];
            double In;

            if ( fmod(Time,TAU) < TINF ){
            In =  DOSE/TINF;
            }
            else {
            In  =  0;
            }

            dA[0] = In - (CL/V)*A[0];
            return List::create(dA);
            }')

PK_1comp_ode_inf_rcpp_ff <- function(model_switch, xt, parameters, poped.db){
  with(as.list(parameters),{
    A_ini <- c(A1=0)

    times_xt <- drop(xt)
    times <- c(0,times_xt) ## add extra time for start of integration
    times <- sort(times) 
    times <- unique(times) # remove duplicates

    out <- ode(A_ini, times, one_comp_ode, parameters)#, events = list(data = eventdat))#atol=1e-13,rtol=1e-13)
    y = out[, "A1"]/(V)
    y=y[match(times_xt,out[,"time"])]
    y=cbind(y)
    return(list(y=y,poped.db=poped.db))
  })
}

poped.db.rcpp <- create.poped.database(ff_file = "PK_1comp_ode_inf_rcpp_ff",
                                  fError_file="feps.add.prop",
                                  fg_file="sfg",
                                  groupsize=20,
                                  m=1,      #number of groups
                                  sigma=c(0.04,0.01),
                                  notfixed_sigma = c(1,0),
                                  bpop=c(CL=20,V=72.8), 
                                  d=c(CL=0.09,V=0.09), 
                                  xt=c(1,2,8,12,24),
                                  minxt=0,
                                  maxxt=100,
                                  a=c(DOSE=100,TINF=10,TAU=24))

plot_model_prediction(poped.db.rcpp)
tic(); evaluate_design(poped.db.rcpp); toc()

output <- poped_optim(poped.db.rcpp,opt_xt = T,method=c("LS"),parallel=T)
plot_model_prediction(output$poped.db)

Example 10 on the GitHub repository has a compiled version of IV infusion as well: PopED/inst/examples/ex.10.PKPD.HCV.compiled.R. This example is available on your poped installation at:

system.file("examples", package="PopED")

[Sibojang9]

That code works greatly. I didn't know the PKPD.HCV.compiled example, since I searched the content of the package with the "infusion" keyword but got no hits. Nice package and nice examples !!! Thank you.

[cohanlon2]

Hi Andrew,

I am having trouble implementing a short IV infusion. The rationale for wanting to do this being that an IV ‘bolus’ is actually given over a period of time e.g. 3 minutes. Using the code you posted on 11 Jul 2017, with low values of TINF = 0.5, the model prediction graph seems to miss some of the subsequent doses. If TINF = 1, there doesn’t seem to be this problem. The issue seems to depend on the DOSE/TINF value (and therefore the amount “In”).

Do you have any suggestions getting this short IV infusion to work?

Best, Conor

[andrewhooker]

Hi Conor,

The problem you are having is that the ODE solver in deSolve is taking too large steps and missing the infusions. To correct for this I recommend adding some extra times in the ODE solution when dosing events occur. I also recommend using the hmax option to the ode function to set a maximum step size for the solver. In the code below I also add extra precision to the solution with the atol and rtol options. Lastly, I add some extra points to the plot_model_prediction function so that you see the infusions.

library(PopED)

library(deSolve)

PK_1comp_ode_inf <- function(Time, State, Pars){
  with(as.list(c(State, Pars)), {
    if( Time %% TAU < TINF ){ 
      In  <-  DOSE/TINF
    } else {
      In  <-  0
    }
    dA1 <- In - CL/V*A1
    return(list(c(dA1)))
  })
}

PK_1comp_ode_inf_ff <- function(model_switch, xt, parameters, poped.db){
  with(as.list(parameters),{
    A_ini <- c(A1=0)

    times_xt <- drop(xt)
    event_times <-  seq(from=0,to=max(times_xt),by=TAU)
    times <- c(0,times_xt,event_times) ## add extra time for start of integration
    times <- sort(times) 
    times <- unique(times) # remove duplicates

    out <- ode(A_ini, times, PK_1comp_ode_inf, parameters,atol=1e-8,rtol=1e-8,hmax=0.1)
    y = out[, "A1"]/(V)
    y=y[match(times_xt,out[,"time"])]
    y=cbind(y)
    return(list(y=y,poped.db=poped.db))
  })
}

## -- parameter definition function 
## -- names match parameters in function ff
sfg <- function(x,a,bpop,b,bocc){
  parameters=c( CL=bpop[1]*exp(b[1]),
                V=bpop[2]*exp(b[2]),
                DOSE=a[1],
                TINF=a[2],
                TAU=a[3])
  return( parameters ) 
}

poped.db <- create.poped.database(ff_file = "PK_1comp_ode_inf_ff",
                                  fError_file="feps.add.prop",
                                  fg_file="sfg",
                                  groupsize=20,
                                  m=1,      #number of groups
                                  sigma=c(0.04,0.01),
                                  notfixed_sigma = c(1,0),
                                  bpop=c(CL=20,V=72.8), 
                                  d=c(CL=0.09,V=0.09), 
                                  xt=c(1,2,8,12,24),
                                  minxt=0,
                                  maxxt=100,
                                  a=c(DOSE=100,TINF=0.5,TAU=24))

plot_model_prediction(poped.db,model_num_points = 500)

# slow evaluation time
tic(); evaluate_design(poped.db); toc()

[cohanlon2]

Great thanks for your help!

Conor

[Anks2030]

Hi Andrew, This is very useful. Can you please share your experience of adding

1. a combined dosing design like having a loading IV infusion dose followed by an IV infusion maintenance dose? This will be really helpful
2. a WT as a covariate on CL and V, how can I embed the WT distribution if I have the mean and SD of WT ?
3. Finally, if I have some endogenous baseline Protein levels along with their variability, how and where can I add this BL values on the top of predicted y Thanks in advance, -Ankur