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

andreasmeid · 2023-07-17T09:59:04Z

andreasmeid
Jul 17, 2023

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 · 2023-07-18T03:01:36Z

mattfidler
Jul 18, 2023
Maintainer

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, time=1.5, 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}

0 replies

mattfidler · 2023-07-18T03:05:22Z

mattfidler
Jul 18, 2023
Maintainer

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.

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}

0 replies

mattfidler · 2023-07-18T03:07:34Z

mattfidler
Jul 18, 2023
Maintainer

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

4 replies

andreasmeid Jul 18, 2023
Author

Dear Matt,

wonderful, that hits the spot. I had not suspected a bug either, but rather a user hint. This is well implemented with the infusion on/off by the EVID=5 indicator. And the simulated curve matches the (empirical) observations very well (and identical to the Matlab simulation, which was not shown here). All in all, very pleasing.

Many thanks and best regards

Andreas

mattfidler Jul 18, 2023
Maintainer

Great Andreas!

One note, for infusions you typically use the rate, so my earlier simulations really should have been divided by the duration of infusion; if you do that it doesn't match which means something is a bit wrong with the model? (I made the mistake myself yesterday)

library(tidyverse)
library(dplyr)
library(rxode2)
#> rxode2 2.0.13.9000 using 8 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 (Ubuntu 11.3.0-1ubuntu1~22.04.1) 11.3.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)/1.5 ,#nmol/hr
     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-18 with reprex v2.0.2}

andreasmeid Jul 21, 2023
Author

Dear Matt,

good point, which makes totally sense. With the workaround to switch the infusion on/off with EVID=5, it would also be possible to simulate a bolus administration, where a lot is administered in a short time. Probably the model or perhaps also the Matlab simulation must be checked, because it actually fit exactly before...

Best wishes

Andreas

mattfidler Jul 24, 2023
Maintainer

Hi @andreasmeid

Indeed it is possible, you simply change the dose to a bolus (evid=1). You could even make it first order if you wish (your first example is technically first order with ka=1).

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

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

{{title}}

Replies: 3 comments 4 replies

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{editor}}'s edit

{{editor}}'s edit

{{title}}

{{title}}

{{title}}

{{title}}

{{title}}

Select a reply

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

andreasmeid Jul 17, 2023

Replies: 3 comments · 4 replies

mattfidler Jul 18, 2023 Maintainer

mattfidler Jul 18, 2023 Maintainer

mattfidler Jul 18, 2023 Maintainer

andreasmeid Jul 18, 2023 Author

mattfidler Jul 18, 2023 Maintainer

andreasmeid Jul 21, 2023 Author

mattfidler Jul 24, 2023 Maintainer

andreasmeid
Jul 17, 2023

Replies: 3 comments 4 replies

mattfidler
Jul 18, 2023
Maintainer

mattfidler
Jul 18, 2023
Maintainer

mattfidler
Jul 18, 2023
Maintainer

andreasmeid Jul 18, 2023
Author

mattfidler Jul 18, 2023
Maintainer

andreasmeid Jul 21, 2023
Author

mattfidler Jul 24, 2023
Maintainer