Infusion into minimal PBPK mit INPUT-compartment and rate/dur

[andreasmeid]

Dear Matt,

I would like to give an infusion of 1.5 hours into an ODE system of a minimal PBPK model. Here, the units of the compartments are obviously concentrations (nM) (not masses alone), so that the mg administration would also have to be related to the compartment volume. I therefore thought of an additional INPUT compartment (similar to "depot"), for which I then set the infusion duration of 1.5 hours with the et table with "dur". Unfortunately, this gives me somewhat too low simulation values compared to the Matlab template (see R file attached). Am I making a general error of thinking with the infusion INPUT compartment?

Many thanks and best regards

Andreas infusion_INPUT_rate.zip

[mattfidler]

Ai @andreasmeid

If you want an infusion the form is not quite the same, the form is simply:

d/dt(Lp)= INPUT/Vp  +  ...

With the following ODE:

d/dt(INPUT) = -INPUT

This is no longer a zero order infusion process. It should not depend on the amount in the INPUT compartment.

In rxode2 the easiest way to have an infusion into the compartment is to simply apply the infusion to the compartment itself.

If you want to use the INPUT an an infusion you could do the following:

d/dt(INPUT) = 0

And then use the replacement event in rxdoe2 to simulate the behavior of the infusion:

Tokuda_dosing_D01 <- et(amt=1*m*1000000000/(MW*1000) ,#nmol
                                                   ii=0, time=0, cmt="INPUT", evid=5) %>%
et(amt=0, ii=0, cmt="INPUT", evid=5)

Using that I get the following:

library(tidyverse)
library(dplyr)
library(rxode2)
#> rxode2 2.0.13 using 4 threads (see ?getRxThreads)
#>   no cache: create with `rxCreateCache()`
library(scales)
#> 
#> Attaching package: 'scales'
#> The following object is masked from 'package:purrr':
#> 
#>     discard
#> The following object is masked from 'package:readr':
#> 
#>     col_factor

# Parameters
TMDD1_pars <- c(MW=145531.86, # molar mass of Trastuzamab g/mol
                kon=2.52, # binding constant to target site (rate constant)
                koff=1.26, # release constant from target site
                ksyn=0.83, # synthesis rate of receptor
                kdeg=0.1, # degradation rate
                kint=0.01, # internalization rate of receptor-ligand complex
                CLint=0.01751,
                kel=0.004, # elimination rate from plasma
                Q=300, # distribution to/from tissue
                Vc=3, # volume of distribution
                Vp=3)

# Model
mod_TMDD1 <- rxode2({

  ktp=Q/Vp;
  kpt=Q/Vc;

  d/dt(INPUT) = 0;

  d/dt(Lp)= INPUT/Vp -kon*Lp*Rp + koff*LRp - (kel+kpt)*Lp + ktp*Lt; # free ligand (i.e., mAb) plasma
  d/dt(Rp)= ksyn - kon*Lp*Rp + koff*LRp -kdeg*Rp;# free target (i.e., receptor) in plasma
  d/dt(LRp)= kon*Lp*Rp - (kint + koff)*LRp; # mAb- receptor complex in plasma
  d/dt(Lt)= kpt*Lp - ktp*Lt; # 

  Cc=Lp*MW/1000;

})
#> using C compiler: 'gcc.exe (GCC) 12.2.0'

# State intials  
TMDD1_inits <- c(Lp=0, Rp=3.5, LRp=0, Lt=0)
MW=145531.86 # molar mass of Trastuzamab g/mol
m=70;

# dosing
Tokuda_dosing_D01 <- et(timeUnits = "hour") %>% 
  et(amt=1*m*1000000000/(MW*1000) ,#nmol
     ii=0, time=0, cmt="INPUT", evid=5) %>% 
  et(amt=0 ,#nmol
     ii=0, time=1.5, cmt="INPUT", evid=5) %>%
  add.sampling(seq(0,24*25,1))

reference_data_1mgpkg <- data.frame(time=c(0.05,0.08,0.09,0.27,0.50,1.06,2.05,3.04,7.03,9.97,13.99,21.01),
                                    Cc=c(59,2532,6882,18146,15438,12872,9412,8254,3131,960,301,55),
                                    Dose="1 mpk")

TMDD1_pred_D01 <- mod_TMDD1 %>% rxode2::rxSolve(params=TMDD1_pars, events=Tokuda_dosing_D01, inits=TMDD1_inits) %>% mutate(Dose="1 mpk")

TMDD1_pred_D01 %>% as.data.frame() %>% ggplot(., aes(x=as.numeric(time)/24, y=Cc, group=Dose, color=Dose)) + geom_line() +
  geom_point(data=reference_data_1mgpkg, mapping=aes(x=time, y=Cc, group=Dose, color=Dose)) +
  scale_y_continuous("Concentrations (ng/mL)", limits=c(0.01, 1000000), breaks=c(0.01, 1, 100, 10000, 1000000),
                     labels=trans_format('log10',math_format(10^.x)),
                     trans="log10", expand=c(0,0))+
  scale_x_continuous("Time (Days)", limits=c(0,25), breaks=seq(0,25,5), expand=c(0,0))+
  scale_color_manual(values=c("1 mpk"="red",
                              "2 mpk"="blue",
                              "4 mpk"="magenta",
                              "8 mpk"="green"))
#> Warning: Transformation introduced infinite values in continuous y-axis

^{Created on 2023-07-17 with reprex v2.0.2}

[mattfidler]

This seems to match the points you have better (presumably from matlab). It doesn't match exactly though since you are probably using a different underlying ODE solver.

``` r library(tidyverse) library(dplyr) library(rxode2) #> rxode2 2.0.13 using 4 threads (see ?getRxThreads) #> no cache: create with `rxCreateCache()` library(scales) #> #> Attaching package: 'scales' #> The following object is masked from 'package:purrr': #> #> discard #> The following object is masked from 'package:readr': #> #> col_factor # Parameters TMDD1_pars <- c(MW=145531.86, # molar mass of Trastuzamab g/mol kon=2.52, # binding constant to target site (rate constant) koff=1.26, # release constant from target site ksyn=0.83, # synthesis rate of receptor kdeg=0.1, # degradation rate kint=0.01, # internalization rate of receptor-ligand complex CLint=0.01751, kel=0.004, # elimination rate from plasma Q=300, # distribution to/from tissue Vc=3, # volume of distribution Vp=3) # Model mod_TMDD1 <- rxode2({ ktp=Q/Vp; kpt=Q/Vc; d/dt(INPUT) = -INPUT; d/dt(Lp)= INPUT/Vp -kon*Lp*Rp + koff*LRp - (kel+kpt)*Lp + ktp*Lt; # free ligand (i.e., mAb) plasma d/dt(Rp)= ksyn - kon*Lp*Rp + koff*LRp -kdeg*Rp;# free target (i.e., receptor) in plasma d/dt(LRp)= kon*Lp*Rp - (kint + koff)*LRp; # mAb- receptor complex in plasma d/dt(Lt)= kpt*Lp - ktp*Lt; # Cc=Lp*MW/1000; }) # State intials TMDD1_inits <- c(Lp=0, Rp=3.5, LRp=0, Lt=0) MW=145531.86 # molar mass of Trastuzamab g/mol m=70; # dosing Tokuda_dosing_D01 <- et(timeUnits = "hour") %>% et(amt=1*m*1000000000/(MW*1000) ,#nmol ii=0, time=0, cmt="INPUT", dur=1.5) %>% add.sampling(seq(0,24*25,1)) reference_data_1mgpkg <- data.frame(time=c(0.05,0.08,0.09,0.27,0.50,1.06,2.05,3.04,7.03,9.97,13.99,21.01), Cc=c(59,2532,6882,18146,15438,12872,9412,8254,3131,960,301,55), Dose="1 mpk") TMDD1_pred_D01 <- mod_TMDD1 %>% rxode2::rxSolve(params=TMDD1_pars, events=Tokuda_dosing_D01, inits=TMDD1_inits) %>% mutate(Dose="1 mpk") TMDD1_pred_D01 %>% as.data.frame() %>% ggplot(., aes(x=as.numeric(time)/24, y=Cc, group=Dose, color=Dose)) + geom_line() + geom_point(data=reference_data_1mgpkg, mapping=aes(x=time, y=Cc, group=Dose, color=Dose)) + scale_y_continuous("Concentrations (ng/mL)", limits=c(0.01, 1000000), breaks=c(0.01, 1, 100, 10000, 1000000), labels=trans_format('log10',math_format(10^.x)), trans="log10", expand=c(0,0))+ scale_x_continuous("Time (Days)", limits=c(0,25), breaks=seq(0,25,5), expand=c(0,0))+ scale_color_manual(values=c("1 mpk"="red", "2 mpk"="blue", "4 mpk"="magenta", "8 mpk"="green")) #> Warning: Transformation introduced infinite values in continuous y-axis ```  ^{Created on 2023-07-17 with [reprex v2.0.2](https://reprex.tidyverse.org)}

[mattfidler]

I'm moving this to the discussion since I don't believe it is a bug. If I am missing something let me know.