search

metrumresearchgroup / mrgsolve

Simulate from ODE-based population PK/PD and QSP models in R
https://mrgsolve.org
GNU General Public License v2.0
132 stars 36 forks source link

Transit model #1098

Closed EilanAlhersh closed 11 months ago

EilanAlhersh commented 1 year ago

Hello Kyle

we are trying to run this model. the model is built but when i pipe it to my dataset all the cent, peri1 and peri2 and PC are 0!. I'm not sure if the problem is with my dataset or with the model.

here is my model code:

code <- '

$PARAM

ALB = 38.00 AGE = 52.00 WT = 72.90 BIO = 1 F1 = 0

$SET delta=0.1

$CMT @annotated

cent : Central Compartement (mg); PERI1 : Peripheral compartment1 (mg); PERI2 : Peripheral compartement2 (mg);

$GLOBAL int NDOSE = 0; double dosetime[100]; double dose[100];

$MAIN double TVCL = THETA(1); double TVCL1 = CLCOV*TVCL;

double TVV1 = THETA(2); double TVV11 = V1COV*TVV1;

double TVNTR = THETA(3);

double TVQ1 = THETA(4); double TVV2 = THETA(5); double TVQ2 = THETA(6); double TVV3 = THETA(7); double TVMTT = THETA(8);

double CL = TVCL1 exp(ETA(1)); double V1 = TVV11 exp(ETA(2)); double NTR = TVNTR * exp(ETA(4));

double MTT = TVMTT * exp(ETA(3));

double Q1 = TVQ1; double V2 = TVV2; double Q2 = TVQ2; double V3 = TVV3;

double K=CL/V1; double K12 = Q1/V1; double K21 = Q1/V2; double K13 = Q2/V1; double K31 = Q2/V3;

double V1WT = ( 1 + THETA(12)*(WT - 72.90)); double V1COV = V1WT;

double CLWT = ( 1 + THETA(11)(WT - 72.90)); double CLALB = ( 1 + THETA(10)(ALB - 38.00)); double CLAGE = ( 1 + THETA(9)(AGE - 52.00)); double CLCOV = CLAGECLALB*CLWT;

if(NEWIND < 2) NDOSE = 0;

if(self.amt > 0 && self.cmt==1) { NDOSE = NDOSE + 1; dosetime[NDOSE] = self.time; dose[NDOSE] = self.amt; }

double KTR = (NTR+1) / MTT; double NFAC = exp(lgamma(NTR+1)); double KINPT = BIO * pow(KTR,(NTR+1)) / NFAC;

$ODE

double INP = 0; int i = 0; while(i <= NDOSE) { double IPT = 0; if(SOLVERTIME >= dosetime[i]) { double delta = SOLVERTIME - dosetime[i]; IPT = dose[i] pow(delta, NTR) exp(-KTR * delta);
} INP = INP + IPT; ++i; }

dxdt_cent=INP- Kcent-K12cent+ K21PERI1- K13cent + K31PERI2; dxdt_PERI1 = K12cent - K21PERI1; dxdt_PERI2= K13cent - K31*PERI2;

$THETA @annotated

10.4 : CL (L/h) 233 : V1 (L) 2.61 : NTR (1/h) 9.48 : Q1 101 : V2 1.45 : Q2 73.9 : V3 0.847 : MTT -0.00577 : CLAGE1 0.0164 : CLALB1 0.00799 : CLWT1 0.00970 : V1WT1

$OMEGA @block 0.210 0.177 0.320 0 0 0.333 0 0 0 1.03

$SIGMA 0

$TABLE @annotated capture PC = cent/V1; capture Aucinf = 90 / CL;

'

mod <- mcode("transit", code, delta = 0.1) model1 <- mod %>% data_set(mydataset) %>% mrgsim_df()

and here is how my dataset looks like.

[

Capture

](url)

Thank you so much

kylebaron commented 1 year ago

Hi @EilanAlhersh -

Some updates to your model below.

Let me know if this is what you are aiming at.

Kyle

library(mrgsolve)
#> 
#> Attaching package: 'mrgsolve'
#> The following object is masked from 'package:stats':
#> 
#>     filter
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union

code <- '
$PARAM
ALB = 38.00
AGE = 52.00
WT = 72.90
BIO = 1
F1 = 0

$SET delta=0.1

$CMT @annotated
depot  : depot (mg)
cent : Central Compartement (mg);
PERI1 : Peripheral compartment1 (mg);
PERI2 : Peripheral compartement2 (mg);

$GLOBAL
int NDOSE = 0;
double dosetime[100];
double dose[100];

$MAIN
double TVCL = THETA(1);
double TVCL1 = CLCOV*TVCL;

double TVV1 = THETA(2);
double TVV11 = V1COV*TVV1;

double TVNTR = THETA(3);

double TVQ1 = THETA(4);
double TVV2 = THETA(5);
double TVQ2 = THETA(6);
double TVV3 = THETA(7);
double TVMTT = THETA(8);

double CL = TVCL1 * exp(ETA(1));
double V1 = TVV11 * exp(ETA(2));
double NTR = TVNTR * exp(ETA(4));

double MTT = TVMTT * exp(ETA(3));

double Q1 = TVQ1;
double V2 = TVV2;
double Q2 = TVQ2;
double V3 = TVV3;

double K   = CL/V1;
double K12 = Q1/V1;
double K21 = Q1/V2;
double K13 = Q2/V1;
double K31 = Q2/V3;

double V1WT  = (1 + THETA(12)*(WT - 72.90));
double V1COV = V1WT;

double CLWT = ( 1 + THETA(11)*(WT - 72.90));
double CLALB = ( 1 + THETA(10)*(ALB - 38.00));
double CLAGE = ( 1 + THETA(9)*(AGE - 52.00));
double CLCOV = CLAGE*CLALB*CLWT;

if(NEWIND < 2) NDOSE = 0;

if(self.amt > 0 && self.cmt==1) {
NDOSE = NDOSE + 1;
dosetime[NDOSE] = self.time;
dose[NDOSE] = self.amt;
}

double KTR = (NTR+1) / MTT;
double NFAC = exp(lgamma(NTR+1));
double KINPT = BIO * pow(KTR,(NTR+1)) / NFAC;

F_depot = F1;

$ODE

double INP = 0;
int i = 0;
while(i <= NDOSE) {
double IPT = 0;
if(SOLVERTIME >= dosetime[i]) {
double delta = SOLVERTIME - dosetime[i];
IPT = dose[i] * pow(delta, NTR) * exp(-KTR * delta);
}
INP = INP + IPT;
++i;
}

dxdt_cent= KTR * depot - K*cent-K12*cent+ K21*PERI1- K13*cent + K31*PERI2;
dxdt_PERI1 = K12*cent - K21*PERI1;
dxdt_PERI2= K13*cent - K31*PERI2;
dxdt_depot = KINPT * INP - KTR * depot;

$THETA @annotated

10.4 : CL (L/h)
233 : V1 (L)
2.61 : NTR (1/h)
9.48 : Q1
101 : V2
1.45 : Q2
73.9 : V3
0.847 : MTT
-0.00577 : CLAGE1
0.0164 : CLALB1
0.00799 : CLWT1
0.00970 : V1WT1

$OMEGA @block
0.210
0.177 0.320
0 0 0.333
0 0 0 1.03

$SIGMA
0

$TABLE @annotated
capture PC = cent/V1;
capture Aucinf = 90 / CL;

'
data <- expand.ev(amt = 90000, WEIGHT = 82.1, ALB = 42)
mod <- mcode("transit", code, delta = 0.1) %>% zero_re()
#> Building transit ...
#> done.
mod %>% data_set(data) %>% mrgsim() %>% plot()


idata <- data.frame(THETA8 = seq(1, 8))

mod %>% 
  ev(amt = 9000) %>% 
  idata_set(idata) %>% 
  mrgsim(end = 24) %>% plot()


idata <- data.frame(THETA3 = seq(1, 8))

mod %>% 
  ev(amt = 9000) %>% 
  idata_set(idata) %>% 
  mrgsim(end = 4) %>% plot()

Created on 2023-07-24 with reprex v2.0.2

EilanAlhersh commented 1 year ago

thank you so much. But it didn't work. I want to use my own dataset, my aim is to keep all the unique ID and their individualized covariates and PK parameters for further analysis. The model is built successfully, but when I pipe it to my dataset I still get the 0 for all the simulated parameters.

I believe something wrong is related to my dataset, however I can't figure it out. The code of my dataset is

mydataset = study_renal %>% rename(DV = CON, AGE = AGEIC, ID = USUBJID) %>% mutate(amt = ifelse(TIME == 0, 90000, 0), ID = as.numeric(factor(ID, levels = unique(ID))), evid = ifelse(amt > 0, 1, 0 ), cmt = 1)

kylebaron commented 1 year ago

Hi @EilanAlhersh -

If you'd like help with the data set, please share a part of the data set along with the model code you are using in a way that it reproduces the issue.

Kyle

EilanAlhersh commented 1 year ago

Attached is the dataset I'm using. Her is all the code I used in this model:

library(mrgsolve)

code <- '

$PARAM

ALB = 38.00 AGE = 52.00 WT = 72.90 BIO = 1 F1 = 0

$SET delta=0.1

$CMT @annotated

depot : depot (mg); cent : Central Compartement (mg); PERI1 : Peripheral compartment1 (mg); PERI2 : Peripheral compartement2 (mg);

$GLOBAL int NDOSE = 0; double dosetime[100]; double dose[100];

$MAIN double TVCL = THETA(1); double TVCL1 = CLCOV*TVCL;

double TVV1 = THETA(2); double TVV11 = V1COV*TVV1;

double TVNTR = THETA(3);

double TVQ1 = THETA(4); double TVV2 = THETA(5); double TVQ2 = THETA(6); double TVV3 = THETA(7); double TVMTT = THETA(8);

double CL = TVCL1 exp(ETA(1)); double V1 = TVV11 exp(ETA(2)); double NTR = TVNTR * exp(ETA(4));

double MTT = TVMTT * exp(ETA(3));

double Q1 = TVQ1; double V2 = TVV2; double Q2 = TVQ2; double V3 = TVV3;

double K=CL/V1; double K12 = Q1/V1; double K21 = Q1/V2; double K13 = Q2/V1; double K31 = Q2/V3;

double V1WT = ( 1 + THETA(12)*(WT - 72.90)); double V1COV = V1WT;

double CLWT = ( 1 + THETA(11)(WT - 72.90)); double CLALB = ( 1 + THETA(10)(ALB - 38.00)); double CLAGE = ( 1 + THETA(9)(AGE - 52.00)); double CLCOV = CLAGECLALB*CLWT;

if(NEWIND < 2) NDOSE = 0;

if(self.amt > 0 && self.cmt==1) { NDOSE = NDOSE + 1; dosetime[NDOSE] = self.time; dose[NDOSE] = self.amt; }

double KTR = (NTR+1) / MTT; double NFAC = exp(lgamma(NTR+1)); double KINPT = BIO * pow(KTR,(NTR+1)) / NFAC;

F_depot = F1;

$ODE

double INP = 0; int i = 0; while(i <= NDOSE) { double IPT = 0; if(SOLVERTIME >= dosetime[i]) { double delta = SOLVERTIME - dosetime[i]; IPT = dose[i] pow(delta, NTR) exp(-KTR * delta);
} INP = INP + IPT; ++i; }

dxdt_cent= KTR depot - Kcent-K12cent+ K21PERI1- K13cent + K31PERI2; dxdt_PERI1 = K12cent - K21PERI1; dxdt_PERI2= K13cent - K31PERI2; dxdt_depot = KINPT INP - KTR depot;

$THETA @annotated

10.4 : CL (L/h) 233 : V1 (L) 2.61 : NTR (1/h) 9.48 : Q1 101 : V2 1.45 : Q2 73.9 : V3 0.847 : MTT -0.00577 : CLAGE1 0.0164 : CLALB1 0.00799 : CLWT1 0.00970 : V1WT1

$OMEGA @block 0.210 0.177 0.320 0 0 0.333 0 0 0 1.03

$SIGMA 0

$TABLE @annotated capture PC = cent/V1; capture Aucinf = 90 / CL;

'

mydataset = study_renal %>% rename(DV = CON, AGE = AGEIC, ID = USUBJID) %>% mutate(amt = ifelse(TIME == 0, 90000, 0), ID = as.numeric(factor(ID, levels = unique(ID))), evid = ifelse(amt > 0, 1, 0 ), cmt = 1)

mod <- mcode("transit", code, delta = 0.1) %>% zero_re() model1 <- mod %>% data_set(mydataset) %>% mrgsim_df()

Thank you so much

example.csv

kylebaron commented 1 year ago

Let's try that again; here's the simulation from the data you posted. 

library(mrgsolve)
#> 
#> Attaching package: 'mrgsolve'
#> The following object is masked from 'package:stats':
#> 
#>     filter
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union

code <- '
$PARAM
ALB = 38.00
AGE = 52.00
WT = 72.90
BIO = 1
F1 = 0

$SET delta=0.1

$CMT @annotated
depot  : depot (mg)
cent : Central Compartement (mg);
PERI1 : Peripheral compartment1 (mg);
PERI2 : Peripheral compartement2 (mg);

$GLOBAL
int NDOSE = 0;
double dosetime[100];
double dose[100];

$MAIN
double TVCL = THETA(1);
double TVCL1 = CLCOV*TVCL;

double TVV1 = THETA(2);
double TVV11 = V1COV*TVV1;

double TVNTR = THETA(3);

double TVQ1 = THETA(4);
double TVV2 = THETA(5);
double TVQ2 = THETA(6);
double TVV3 = THETA(7);
double TVMTT = THETA(8);

double CL = TVCL1 * exp(ETA(1));
double V1 = TVV11 * exp(ETA(2));
double NTR = TVNTR * exp(ETA(4));

double MTT = TVMTT * exp(ETA(3));

double Q1 = TVQ1;
double V2 = TVV2;
double Q2 = TVQ2;
double V3 = TVV3;

double K   = CL/V1;
double K12 = Q1/V1;
double K21 = Q1/V2;
double K13 = Q2/V1;
double K31 = Q2/V3;

double V1WT  = (1 + THETA(12)*(WT - 72.90));
double V1COV = V1WT;

double CLWT = ( 1 + THETA(11)*(WT - 72.90));
double CLALB = ( 1 + THETA(10)*(ALB - 38.00));
double CLAGE = ( 1 + THETA(9)*(AGE - 52.00));
double CLCOV = CLAGE*CLALB*CLWT;

if(NEWIND < 2) NDOSE = 0;

if(self.amt > 0 && self.cmt==1) {
NDOSE = NDOSE + 1;
dosetime[NDOSE] = self.time;
dose[NDOSE] = self.amt;
}

double KTR = (NTR+1) / MTT;
double NFAC = exp(lgamma(NTR+1));
double KINPT = BIO * pow(KTR,(NTR+1)) / NFAC;

F_depot = F1;

$ODE

double INP = 0;
int i = 0;
while(i <= NDOSE) {
double IPT = 0;
if(SOLVERTIME >= dosetime[i]) {
double delta = SOLVERTIME - dosetime[i];
IPT = dose[i] * pow(delta, NTR) * exp(-KTR * delta);
}
INP = INP + IPT;
++i;
}

dxdt_cent= KTR * depot - K*cent-K12*cent+ K21*PERI1- K13*cent + K31*PERI2;
dxdt_PERI1 = K12*cent - K21*PERI1;
dxdt_PERI2= K13*cent - K31*PERI2;
dxdt_depot = KINPT * INP - KTR * depot;

$THETA @annotated

10.4 : CL (L/h)
233 : V1 (L)
2.61 : NTR (1/h)
9.48 : Q1
101 : V2
1.45 : Q2
73.9 : V3
0.847 : MTT
-0.00577 : CLAGE1
0.0164 : CLALB1
0.00799 : CLWT1
0.00970 : V1WT1

$OMEGA @block
0.210
0.177 0.320
0 0 0.333
0 0 0 1.03

$SIGMA
0

$TABLE @annotated
capture PC = cent/V1;
capture Aucinf = 90 / CL;

'
data <- expand.ev(amt = 90000, WEIGHT = 82.1, ALB = 42)
mod <- mcode("transit", code, delta = 0.1) %>% zero_re()
#> Building transit ...
#> done.
mod %>% data_set(data) %>% mrgsim() %>% plot()


# idata <- data.frame(THETA8 = seq(1, 8))
# 
# mod %>% 
#   ev(amt = 9000) %>% 
#   idata_set(idata) %>% 
#   mrgsim(end = 24) %>% plot()
# 
# 
# idata <- data.frame(THETA3 = seq(1, 8))
# 
# mod %>% 
#   ev(amt = 9000) %>% 
#   idata_set(idata) %>% 
#   mrgsim(end = 4) %>% plot()

You need to make the data names consistent

See ?uctran

data <- data.table::fread("~/helpr/example.csv", na.strings = ".")
data <- uctran(data)
mrgsim(mod, data) %>% filter_sims(TIME <= 4) %>% plot()

packageVersion("mrgsolve")
#> [1] '1.0.9.9000'

Created on 2023-07-25 with reprex v2.0.2

EilanAlhersh commented 1 year ago

It works! OMG, I just rename TIME to time and it works. how sensitive it is.

Thank you so much

kylebaron commented 1 year ago

There is a check of incoming data to ensure consistency of certain data names such as amt / AMT or cmt / CMT and it will error if there is ambiguity. TIME should be in that check, but it wasn't. So it's a bug that will get fixed.