Adding covariates into model

[meensrinivasan]

Hello,

There are a couple of questions I had while coding a model using RxODE.

1.) How to account for weight/BSA based dosing in add.dosing? eg, if I want to simulate 20 mg/m2 dose IV infusion. Now since this depends on BSA, is there any way I need not hard code it and say, get the input from the user? Say in a shiny based app?

ev$add.dosing(dose = 20, rate = 120, nbr.doses = 1)

Here, ideally the rate would depend on the dose, so is there any way to specify that?

2.) My model has continuous (eg, weight) and categorical (eg, sex) covariates, I would like to get this input via the user too. How can we incorporate covariates into the model? The example on the RxODE cran page does mention the input for a time varying covariate, but I was wondering if there is an example I could refer to?

my example-

say I would like to code clearance as a function of Creat clearance, weight and gender

CL=a(WT/70)^0.75(CRCL)^b*c^SEX

a,b and c being theta's for population CL, CRCL on CL and sex on clearance- a,b and c are known. How would I code this?

I appreciate any help on this matter!

Thank you :)

[mattfidler]

Note for these questions I will use the latest stable, which you can always get from devtools:

devtools::install_github("nlmixrdevelopment/RxODE",ref="v0.6-2")

For the first issue:

How to account for weight/BSA based dosing in add.dosing?

You can use standard R to normalize weight or dose into rate:

# Dose is weight normalized
ev$add.dosing(dose = 20/wt, rate = 120, nbr.doses = 1)
# Specify duration instead via the rate construct
dose <- 20;
ev$add.dosing(dose=dose,rate=dose/4);

How can we incorporate covariates into the model?

Take the example model:

mod <- RxODE({
   C2 = centr/V2;
   C3 = peri/V3;
   d/dt(depot) =-KA*depot;
   d/dt(centr) = KA*depot - CL*C2 - Q*C2 + Q*C3;
   d/dt(peri)  =                    Q*C2 - Q*C3;
   d/dt(eff)  = Kin - Kout*(1-C2/(EC50+C2))*eff;
   eff(0) = 1
})

For me this gives the following model:

> mod
RxODE model named "rx_ba412273210f8e8f956119e1940f1949" (ready to run).
States: depot, centr, peri, eff
Params: V2, V3, KA, CL, Q, Kin, Kout, EC50

Note that CL is a parameter in the model.

Say you wanted to add CL=a*(WT/70)^0.75*(CRCL)^b*c^SEX you could do so easily by the following:

mod2 <- RxODE({
   CL=a*(WT/70)^0.75*(CRCL)^b*c^SEX
   C2 = centr/V2;
   C3 = peri/V3;
   d/dt(depot) =-KA*depot;
   d/dt(centr) = KA*depot - CL*C2 - Q*C2 + Q*C3;
   d/dt(peri)  =                    Q*C2 - Q*C3;
   d/dt(eff)  = Kin - Kout*(1-C2/(EC50+C2))*eff;
   eff(0) = 1
})

Note the new model no longer has CL but has the parameters you may wish to input a, WT, etc.

> mod2
RxODE model named "rx_e1bb377ca885955da251f172827fe4d8" (ready to run).
States: depot, centr, peri, eff
Params: a, WT, CRCL, b, c, SEX, V2, V3, KA, Q, Kin, Kout, EC50

If you type summary(mod2) you get more model information:

> summary(mod)
RxODE model named "rx_ba412273210f8e8f956119e1940f1949" (ready to run).
DLL: c:/.../rx_ba412273210f8e8f956119e1940f1949_x64.dll

User supplied parameters:
  V2   V3   KA   CL    Q  Kin Kout EC50 
  NA   NA   NA   NA   NA   NA   NA   NA 

User initial conditions:
depot centr  peri   eff 
    0     0     0     0 

Compartments:
  cmt=1   cmt=2   cmt=3   cmt=4 
"depot" "centr"  "peri"   "eff" 

Calculated Variables:
[1] "C2" "C3"

Model:
    C2 = centr/V2
    C3 = peri/V3
    d/dt(depot) = -KA * depot
    d/dt(centr) = KA * depot - CL * C2 - Q * C2 + Q * C3
    d/dt(peri) = Q * C2 - Q * C3
    d/dt(eff) = Kin - Kout * (1 - C2/(EC50 + C2)) * eff

If these are all known, you can solve them as you would any other RxODE model, by first setting up an event table and then solving the subject:

et <- eventTable();
et$add.sampling(seq(0, 7, length.out=200));
et$add.dosing(20, start.time=0);

theta <- 
   c(KA=2.94E-01, a=1.86E+01, V2=4.02E+01, # central 
     Q=1.05E+01,  V3=2.97E+02,              # peripheral
     Kin=1, Kout=1, EC50=200,               # effects
    WT = 80, CRCL=120,  b=0.125,  c=1.1, SEX=1) # Covariates

x <- solve(mod2, theta, et);

This gives

> x
Solved RxODE object
Dll: c:/SVN/Wenping/RxODE/R/rx_e1bb377ca885955da251f172827fe4d8_x64.dll

Parameters:
      a      WT    CRCL       b       c     SEX      V2      V3      KA       Q 
 18.600  80.000 120.000   0.125   1.100   1.000  40.200 297.000   0.294  10.500 
    Kin    Kout    EC50 
  1.000   1.000 200.000 

Initial Conditions:
depot centr  peri   eff 
    0     0     0     1 

First part of data:
# A tibble: 200 x 8
    time depot centr     peri   eff    CL      C2         C3
   <dbl> <dbl> <dbl>    <dbl> <dbl> <dbl>   <dbl>      <dbl>
1 0       20.0 0     0         1.00  41.1 0       0         
2 0.0352  19.8 0.201 0.000932  1.00  41.1 0.00500 0.00000314
3 0.0704  19.6 0.391 0.00366   1.00  41.1 0.00974 0.0000123 
4 0.106   19.4 0.571 0.00808   1.00  41.1 0.0142  0.0000272 
5 0.141   19.2 0.741 0.0141    1.00  41.1 0.0184  0.0000475 
6 0.176   19.0 0.901 0.0216    1.00  41.1 0.0224  0.0000728 
# ... with 194 more rows

Plotting can give

library(ggplot2)
ggplot(x,aes(time,eff))+geom_line()

And

ggplot(x,aes(time,centr))+geom_line()

If however, one of your covariates are time-varying you can also solve. For instance if CRCL is a time-varying covariate in your dataset. Then provide that as a dataest to the covs argument when solving:

crcl.new <- data.frame(CRCL=seq(88, 128, length.out=200));

theta <- theta[names(theta) != "CRCL"]
x <- solve(mod2, theta, et,covs=crcl.new);

This gives:

> x
Solved RxODE object
Dll: c:/SVN/Wenping/RxODE/R/rx_e1bb377ca885955da251f172827fe4d8_x64.dll

Parameters:
      a      WT       b       c     SEX      V2      V3      KA       Q     Kin 
 18.600  80.000   0.125   1.100   1.000  40.200 297.000   0.294  10.500   1.000 
   Kout    EC50 
  1.000 200.000 

Time Varying Covariates
CRCL

Initial Conditions:
depot centr  peri   eff 
    0     0     0     1 

First part of data:
# A tibble: 200 x 8
    time depot centr     peri   eff    CL      C2         C3
   <dbl> <dbl> <dbl>    <dbl> <dbl> <dbl>   <dbl>      <dbl>
1 0       20.0 0     0         1.00  39.6 0       0         
2 0.0352  19.8 0.201 0.000933  1.00  39.6 0.00501 0.00000314
3 0.0704  19.6 0.392 0.00366   1.00  39.6 0.00975 0.0000123 
4 0.106   19.4 0.572 0.00809   1.00  39.6 0.0142  0.0000273 
5 0.141   19.2 0.743 0.0141    1.00  39.6 0.0185  0.0000476 
6 0.176   19.0 0.904 0.0217    1.00  39.6 0.0225  0.0000730 
# ... with 194 more rows
> 

The covariates are interpolated by linear interpolation when needed by default. NONMEM uses the next observation caried backward (NOCB) and Monolix uses last observation carried forward (LOCF)

This can be solved by the following options:

x <- solve(mod2, theta, et,covs=crcl.new, covs_interpolation="nocb"); # NONMEM style
x <- solve(mod2, theta, et,covs=crcl.new, covs_interpolation="locf"); # Monolix style

[mattfidler]

The development version is working on solving all of this based on a single dataest, like NONMEM (FYI)

[meensrinivasan]

Thank you so much @mattfidler this explanation was super clear and helpful. Just tried it!! 🥇

1.) I tried my entire model, which has covariates on CL, V1, V2 and Q -specified in the mod parameter like you mentioned and specifying the theta values, I just need to find a way to not hard code it (like in the params/theta variable) and rather use the user-defined values, say from a shiny server.

2.) My simulated concentrations are around only 1/4th the values that the original dataset contains, any idea what I could have missed? I've checked the units.

Thank you so much for the detailed explanation!! ^_^ Appreciate your help!

[mattfidler]

I just need to find a way to not hard code it

Functions are the way:

myFun <- function(SEX){
et <- eventTable();
et$add.sampling(seq(0, 7, length.out=200));
et$add.dosing(20, start.time=0);

theta <- 
   c(KA=2.94E-01, a=1.86E+01, V2=4.02E+01, # central 
     Q=1.05E+01,  V3=2.97E+02,              # peripheral
     Kin=1, Kout=1, EC50=200,               # effects
    WT = 80, CRCL=120,  b=0.125,  c=1.1, SEX=SEX) # Covariates

x <- solve(mod2, theta, et);
return(x)
}

My simulated concentrations are around only 1/4th the values that the original dataset contains

Often units or problems in dosing are to blame. I would make sure your rate equations make sense.

[meensrinivasan]

Thanks @mattfidler for the example.

I tried the same in my shiny server code. and I got the error that "Object x not found"

Here is my code

#Input options
DOSE <- input$CDOSE
SEX <- input$ISM
CRCL <- input$CRCL
ALB <- input$ALB
WT <- input$Weight

#Differential equations describing model
DES <- function(WT,ISM,ALB, CRCL) {

  ev <- eventTable();
  for (i in 1:n)ev$add.dosing(dose = DOSE, rate = DOSE*6, nbr.doses = 1);
  for (i in 1:n)ev$add.sampling(0:24);

  params <- c(thCL = 5, thCRCLonCL=0.5, thISMonCL=1.5,
              thV1 = 10, thALBonV1=-2, 
              thQ = 5, 
              thV2 = 10, 
              ALB=ALB, CRCL=CRCL, ISM=ISM, WT=WT)            

  ode <- "
  CL<-thCL*(CRCL/100)^thCRCLonCL*thISMonCL^ISM
  V1<-thV1*(WT/70)*(ALB/3)^thALBonV1
  V2<-thV2*(WT/70)
  Q<-thQ*(WT/70)^0.75
  C1 = A1/V1;   
  C2 = A2/V2;
  d/dt(A1) = -(Q/V1)*A1-(CL/V1)*A1+(Q/V2)*A2;
  d/dt(A2) = (Q/V1)*A1 -(Q/V2)*A2;
  "
  mod1 <- RxODE(model = ode, modName = "mod1")
  x <- solve(mod1, params, ev)
  return(x)}

})

output$plotCONC <- renderPlot({
  ggplot(x,aes(time,C1))+geom_line()
})
output$DOSE <- renderText({
  paste("Dose to be administered =", signif(input$BSA*500, digits = 3) ,"mg")

})  

output$RATE <- renderText({
  paste("Infusion rate =", signif(input$CDOSE*6, digits = 3) ,"mg/hr")

})

})

Anything that I've missed? Thank you :)

[mattfidler]

My guess you are facing an R scoping issue; In my opinion, the DES function should run outside shiny and then integrated afterward:

DES <- function(WT,ISM,ALB, CRCL,n=4,DOSE=200) {

  ev <- eventTable();
  for (i in 1:n)ev$add.dosing(dose = DOSE, rate = DOSE*6, nbr.doses = 1, start.time=0);
  for (i in 1:n)ev$add.sampling(0:24);

  params <- c(thCL = 5, thCRCLonCL=0.5, thISMonCL=1.5,
              thV1 = 10, thALBonV1=-2, 
              thQ = 5, 
              thV2 = 10, 
              ALB=ALB, CRCL=CRCL, ISM=ISM, WT=WT)            

  ode <- "
  CL<-thCL*(CRCL/100)^thCRCLonCL*thISMonCL^ISM
  V1<-thV1*(WT/70)*(ALB/3)^thALBonV1
  V2<-thV2*(WT/70)
  Q<-thQ*(WT/70)^0.75
  C1 = A1/V1;   
  C2 = A2/V2;
  d/dt(A1) = -(Q/V1)*A1-(CL/V1)*A1+(Q/V2)*A2;
  d/dt(A2) = (Q/V1)*A1 -(Q/V2)*A2;
  "
  mod1 <- RxODE(model = ode, modName = "mod1")
  x <- solve(mod1, params, ev)
  return(x)}

Will return the solved parameters. However, this may be a shiny issue.

[mattfidler]

I'm a little unclear the purpose of the n variable, by the way. The events seem to be piling on top of another.

[mattfidler]

It may also be better to create the model in your shiny model setup and then solve it interactively.

[meensrinivasan]

I kept the n variable so that one can select the number of individuals to conduct the simulation on through the ui.r Thank you so much for your help !

[mattfidler]

So you got it working? Congratulations

[meensrinivasan]

Still working on it :)

[mattfidler]

Good Luck :)