Combined PK/PD Model: Definition of error model

[marianklose]

Hi there, thanks for your great work! :) Just wanted to know how I could define a separate error model for a combined PK/PD model. Basically I have defined one comprehensive model file for both PK and PD (neutrophil count).

For the error model I am defining:

error_model = function(f,sigma){
g <- sigma[1] + sigma[2]*f return(g)

for a additive+proportional error. The documentation says that f is the predicted concentration in this function. But how is f defined in case we have two predicted measures of interest? And what if I want to have for example a different error model for PK than for PD? Thanks for your help!

[levenc]

Hi @klosinho00, thank you! Your use case is quite justified, I see what you are trying to do. It is not yet possible to do this in posologyr v1.1.0 yet, but I'll branch tomorrow to see what changes are needed and do some testing. Once this work is done, would you like to help test it before it is merged ?

[marianklose]

Hi, thanks for your quick response and sorry for my late answer. Sure, I am happy to help where I can! Just let me know how we should proceed when you are done.

[levenc]

Dear @klosinho00 , you can now define a separate error model for multi-endpoints models. I hope this works adequately for your PKPD model.

The comparison of posologyr's performance with the individual nlmixr2 estimates is quite satisfactory. If you happen to have other tools that can estimate individual profiles (NONMEM, Monolix...) I would be interested to know how posologyr compares to these references.

First of all, for the time being, you will need to install the version from the multioutput branch:

remotes:: install_github ("posologyr/multioutput")

To get estimations for your integrated PK and PD model, you have to be consistent in your naming convention across:

• your dataset (using the column DVID)
• your residual error models (in a named list)
• the standard deviation of your residual error models (in a named list)

You can define as many residual error models as you like. Each model defined in the named list error_models must have its counterpart in the named list sigma, and the names must be the same as the names defined in the DVID column of the data set.

See the following example, using the warfarin data and model (provided by Tomoo Funaki and Nick Holford) from the nlmixr documentation (https://nlmixr2.org/articles/multiple-endpoints.html).

warfarin PKPD model

mod_warfarin_nlmixr <- list(
  ppk_model   = rxode2::rxode({
    ktr <- exp(THETA_ktr + ETA_ktr)
    ka <- exp(THETA_ka + ETA_ka)
    cl <- exp(THETA_cl + ETA_cl)
    v <- exp(THETA_v + ETA_v)
    emax = expit(THETA_emax + ETA_emax)
    ec50 =  exp(THETA_ec50 + ETA_ec50)
    kout = exp(THETA_kout + ETA_kout)
    e0 = exp(THETA_e0 + ETA_e0)
    ##
    DCP = center/v
    PD=1-emax*DCP/(ec50+DCP)
    ##
    effect(0) = e0
    kin = e0*kout
    ##
    d/dt(depot) = -ktr * depot
    d/dt(gut) =  ktr * depot -ka * gut
    d/dt(center) =  ka * gut - cl / v * center
    d/dt(effect) = kin*PD -kout*effect
    ##
    cp = center / v
    pca = effect
  }),
  error_model = list(
    cp = function(f,sigma){
      g <- sigma[1]^2 + (sigma[2]^2)*(f^2)
      return(sqrt(g))
    },
    pca = function(f,sigma){
      g <- sigma[1]^2 + (sigma[2]^2)*(f^2)
      return(sqrt(g))
    }
  ),
  theta = c(THETA_ktr=0.106,
            THETA_ka=-0.087,
            THETA_cl=-2.03,
            THETA_v=2.07,
            THETA_emax=3.4,
            THETA_ec50=0.00724,
            THETA_kout=-2.9,
            THETA_e0=4.57),
  omega = lotri::lotri({ETA_ktr + ETA_ka + ETA_cl + ETA_v + ETA_emax + ETA_ec50 +
      ETA_kout + ETA_e0 ~
      c(1.024695,
        0.00     ,   0.9518403 ,
        0.00     ,   0.00 ,   0.5300943  ,
        0.00    ,    0.00,   0.00,   0.4785394,
        0.00    ,    0.00,   0.00,   0.00,   0.7134424,
        0.00    ,    0.00,   0.00,   0.00,   0.00,   0.7204165,
        0.00    ,    0.00,   0.00,   0.00,   0.00,   0.00,   0.3563706,
        0.00    ,    0.00,   0.00,   0.00,   0.00,   0.00,   0.00,   0.2660827)}),
  sigma       = list(
    cp=c(additive_a = 0.144, proportional_b = 0.15),
    pca=c(additive_a = 3.91, proportional_b = 0.0)
    )
  )

data: first subject from the warfarin dataset

> warf_01 <- data.frame(ID=1,
+                       TIME=c(0.0,0.5,1.0,2.0,3.0,6.0,9.0,12.0,24.0,24.0,36.0,
+                              36.0,48.0,48.0,72.0,72.0,96.0,120.0,144.0),
+                       DV=c(0.0,0.0,1.9,3.3,6.6,9.1,10.8,8.6,5.6,44.0,4.0,27.0,
+                            2.7,28.0,0.8,31.0,60.0,65.0,71.0),
+                       DVID=c("cp","cp","cp","cp","cp","cp","cp","cp","cp","pca",
+                              "cp","pca","cp","pca","cp","pca","pca","pca","pca"),
+                       EVID=c(1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
+                       AMT=c(100,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0))
> warf_01
   ID  TIME   DV DVID EVID AMT
1   1   0.0  0.0   cp    1 100
2   1   0.5  0.0   cp    0   0
3   1   1.0  1.9   cp    0   0
4   1   2.0  3.3   cp    0   0
5   1   3.0  6.6   cp    0   0
6   1   6.0  9.1   cp    0   0
7   1   9.0 10.8   cp    0   0
8   1  12.0  8.6   cp    0   0
9   1  24.0  5.6   cp    0   0
10  1  24.0 44.0  pca    0   0
11  1  36.0  4.0   cp    0   0
12  1  36.0 27.0  pca    0   0
13  1  48.0  2.7   cp    0   0
14  1  48.0 28.0  pca    0   0
15  1  72.0  0.8   cp    0   0
16  1  72.0 31.0  pca    0   0
17  1  96.0 60.0  pca    0   0
18  1 120.0 65.0  pca    0   0
19  1 144.0 71.0  pca    0   0

posologyr can compute the EBE for the combined PKPD model with poso_estim_map()

> map_warf_01  <- poso_estim_map(warf_01,mod_warfarin_nlmixr)
> map_warf_01
$eta
    ETA_ktr      ETA_ka      ETA_cl       ETA_v    ETA_emax    ETA_ec50    ETA_kout      ETA_e0 
-0.43409547 -0.69456986  0.82337519 -0.02306934  0.06905641  0.14183213 -0.19873468 -0.12379873 

$model
── Solved rxode2 object ──
── Parameters ($params): ──
  THETA_ktr     ETA_ktr    THETA_ka      ETA_ka    THETA_cl      ETA_cl     THETA_v       ETA_v  THETA_emax 
 0.10600000 -0.43409547 -0.08700000 -0.69456986 -2.03000000  0.82337519  2.07000000 -0.02306934  3.40000000 
   ETA_emax  THETA_ec50    ETA_ec50  THETA_kout    ETA_kout    THETA_e0      ETA_e0 
 0.06905641  0.00724000  0.14183213 -2.90000000 -0.19873468  4.57000000 -0.12379873 
── Initial Conditions ($inits): ──
 depot    gut center effect 
     0      0      0      0 
── First part of data (object): ──
# A tibble: 1,451 × 18
   time   ktr    ka    cl     v  emax  ec50   kout    e0    DCP    PD   kin     cp   pca depot   gut center
  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl> <dbl>  <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl>  <dbl>
1   0   0.720 0.458 0.299  7.74 0.970  1.16 0.0451  85.3 0      1      3.85 0       85.3 100    0     0    
2   0.1 0.720 0.458 0.299  7.74 0.970  1.16 0.0451  85.3 0.0204 0.983  3.85 0.0204  85.3  93.1  6.79  0.158
3   0.2 0.720 0.458 0.299  7.74 0.970  1.16 0.0451  85.3 0.0785 0.939  3.85 0.0785  85.3  86.6 12.8   0.608
4   0.3 0.720 0.458 0.299  7.74 0.970  1.16 0.0451  85.3 0.170  0.876  3.85 0.170   85.3  80.6 18.1   1.31 
5   0.4 0.720 0.458 0.299  7.74 0.970  1.16 0.0451  85.3 0.290  0.806  3.85 0.290   85.2  75.0 22.8   2.25 
6   0.5 0.720 0.458 0.299  7.74 0.970  1.16 0.0451  85.3 0.436  0.735  3.85 0.436   85.1  69.8 26.8   3.37 
# ℹ 1,445 more rows
# ℹ 1 more variable: effect <dbl>
# ℹ Use `print(n = ...)` to see more rows

$event
      id  time amt evid
   1:  1   0.0  NA    0
   2:  1   0.0 100    1
   3:  1   0.1  NA    0
   4:  1   0.2  NA    0
   5:  1   0.3  NA    0
  ---                  
1448:  1 144.6  NA    0
1449:  1 144.7  NA    0
1450:  1 144.8  NA    0
1451:  1 144.9  NA    0
1452:  1 145.0  NA    0

You can also plot the observation/time curves for both endpoints

plot(map_warf_01$model,"cp")

plot(map_warf_01$model,"pca")

[marianklose]

Dear @levenc, thanks a lot for working on this! This feature surely adds new helpful functionality to your tool and makes it applicable to a broader range of use cases. I do have access to NONMEM and could compare the MAP estimates of posologyr with NONMEM. It's on my agenda and I will let you know whenever I found time to work on it. Thanks again.

[levenc]

Glad I could help, and thank you for the offer to compare posologyr with NONMEM. I'm curious to see the result if you find the time.