DV missing in output and far too small iiv

andreasmeid commented 3 years ago

Dear Félicien,

sorry for asking again. Now I am trying to apply mapbayr to a published (NONMEM) model. I probably missed something fundamental in the mrgsolve code, because the mapbayest function shows only NA for the variable DV.

Output (11 lines):
  ID time evid cmt amt mdv DV      IPRED       PRED ...
1  1    0    1   1 2.5   1 NA 0.07253891 0.07253863 ...
2  1   12    1   1 2.5   1 NA 0.07550071 0.07549832 ...
3  1   24    1   1 2.5   1 NA 0.07280429 0.07280366 ...
4  1   36    1   1 2.5   1 NA 0.07551154 0.07550914 ...
5  1   48    1   1 2.5   1 NA 0.07280432 0.07280370 ...
6  1   60    1   1 2.5   1 NA 0.07551155 0.07550914 ...

Also, the inter-individual variability (via the omegas) does not seem to work properly; population predictions and individual predictions are practically on top of each other. The estimated omegas must already be used in the code, right (see for example Table 3 of Cirincione B et al. CPT Pharmacometrics Syst Pharmacol 2018;7:728-738)? Or is the are the CV values to be considered for this?

Following, I have attached my code; it would be nice if you could take a critical look at it and maybe spot the errors.

Many thanks and best wishes

Andreas

R Code:

NONMEM_code <- "
$PROB
;----------------------------------------;
;Stage1 Final Model from Cirincione B et al. CPT Pharmacometrics Syst Pharmacol 2018;7:728-738
;----------------------------------------;

$PARAM @annotated
TVKA: 0.471 : Absorption constant morning (1/h) (THETA(1))
TVCL:  3.59 : Clearance (CL/F)
TVV2: 29.8 : Central volume
TVQ   :  1.62 : Intercompartmental clearance
TVV3  : 26.7 : Peripheral volume of distribution

ETA1: 0 : Absorption / random effects
ETA2: 0 : Clearance (L/h)
ETA3: 0 : Central volume (L)
ETA4: 0 : Intercompartmental clearance

$PARAM @annotated @covariates
D_WTB : 70 : Body weight (kg) at baseline
EVENING : 1 : evening dose = 1, morning dose = 0
AF : 1 : NVAF patient indicator (set to 1)
SEX : 0 : 0:Male, 1:Female
AGE : 69 : Age in years at baseline
D_CCRCLB : 70 : Creatinine clearance CCG

$MAIN
double V2 = TVV2*(pow(D_WTB/65,0.806))*(1−0.0217) * exp(ETA3 + ETA(3)) ;
double KA = TVKA * pow((1-0.434), EVENING) * exp(ETA1 + ETA(1)) ; //* pow(0.201, MORNING); // * pow(0.201, MORNING)
double CL = ( (((1.57*(D_CCRCLB/70)) + (2.02* (pow(AGE/69, −0.429))*(1+SEX*(−0.215))) )*(1−0.149)) /TVV2) * exp(ETA2 + ETA(2)) ; // * pow(BW / 70, 1.2) ;
double Q = TVQ * exp(ETA4 + ETA(4)) ;
double V3 = TVV3 ;

// $OMEGA 0.521 0.318 0.169 0.857
$OMEGA @block
0.271 // KA
0.000000 0.101 // CL
0.000000 0.000000 0.0287 // V1
0.000000 0.000000 0.000000 0.867 // Q as mean of omegas from k12 and k21

$SIGMA @labels PROP ADD
0.01 // proportional
0.3 // additive

$CMT @annotated
CENT  : Central compartment (mg/L)[OBS]
GUT  : Dosing compartment  (mg)[ADM]
PERIPH: Peripheral compartment V3 ()

$TABLE
double IPRED = CENT/V2;
double DV = (CENT/V2) *(1+PROP)+ADD;

//while(DV < 0) {
//simeps();
//DV = IPRED*(1+PROP)+ADD;
//}

$PKMODEL ncmt=2, depot=TRUE

$CAPTURE DV CL KA V2 Q V3
"

NONMEN_model <- mcode("NONMEN_model", NONMEM_code)

vorzeitraum <- data.frame(ID=1, time=seq(0,84,12), evid=1, cmt=1, amt=2.5, DV=NA, mdv=1, D_WTB=80, EVENING=rep(c(0,1),8),
                          D_CCRCLB=70,SEX=1, AGE=76)[1:8,]
measure_data <- data.frame(ID = 1, time = c(96,96+12.4,96+19.5), evid = c(1, rep(0,2)), cmt = 1, amt = c(2.5, rep(0,2)),
                      DV = c(NA, 0.19, 0.12), mdv = c(1,0,0), D_WTB  = 80, EVENING=0, D_CCRCLB=70,SEX=1, AGE=76)
NONMEM_data <- rbind(vorzeitraum, measure_data)

test_model %>% data_set(NONMEM_data[1:9,]) %>% mrgsim(start=0,end=120, delta=1)%>% plot(CENT~time|factor(ID))

NONMEM_est <- mapbayest(NONMEN_model, data = NONMEM_data[1:11,])
print(NONMEM_est)
plot(NONMEM_est)

FelicienLL commented 3 years ago

Hi Andreas, My answers so far:

Nothing is wrong to have NA in the DV column of your output, since only administration lines are shown.
In the $OMEGA block, a variance-covariance matrix of inter-individual variability is expected, or, at least its diagonal, i.e. the variances. So, you have to set the values of omega², and not the CV%.
Something feels wrong in the estimation. I suspect a problem with compartments. I strongly advise you to re-order $CMT in your code, with GUT first, then CENTRAL and PERIPH at last. Because in the background I'm not sure how mrgsolve behaves with the $PKMODEL routine when compartment are not properly named... Also, in your data, the drug is both administered and observed in compartment 1 (cmt = 1), which I guess is not correct ?
Once it is done, try to simulate some typical profile to see, at least, if the course of concentration and quantities in each compartment goes as you imagine. From what I see, the drug accumulates a lot... I suspect an error in the formula of CL since just running mrgsim(NONMEN_model) returns a CL of 0.1 (on average), while you specified a typical value of 3.59, so it is maybe unlikely ?

Please let me know if some of these helped - or did not !

Félicien

FelicienLL commented 3 years ago

I also took a look to the published model. Random effects were set on microconstant (K, K12, K21) instead of clearance. For a more accurate approach, I would also suggest you to do the same in your model by defining K, K12 and K21 in $MAIN. Then, just define CL, V3 and Q as function of these pre-defind K,K12 and K2 where IIV was applied. (Or, even easier, implement the microconstants directly in differential equations in $ODE) Félicien

andreasmeid commented 3 years ago

Dear Félicien, thank you very much for these tips;

The actual omegas from the publication are now used
The code was further adapted to the publication, so that now the transfer constants k, k21, and k are included accordingly
By rearranging the compartments (GUT first), cmt=2 had to be set in the dataset
The population predictions now look better (both in mrgsolve and mapbayr); still the iiv variability is literally missing and PPRED and IPRED are practically on top of each other. The etas are vanishingly small.

What can be the reason for this, meaning that practically no iiv exists, although the measured values deviate clearly from the population predictions? I have attached the updated code again and am as always very grateful for every step towards the solution of the problem. Best wishes Andreas

NONMEM_code <- "
$PROB
;----------------------------------------;
;Stage1 Final Model from 
;----------------------------------------;

$PARAM @annotated
TVKA: 0.471 : Absorption constant morning (1/h) 
TVCLR:  1.57 : Renal clearance (L/h)
TVCLNR:  2.02 : Non-renal clearance (L/h)
TVV2: 29.8 : Central volume
TVQ   :  1.62 : Intercompartmental clearance
TVV3  : 26.7 : Peripheral volume of distribution
ETA1: 0 : Absorption / random effects (1/h)
ETA2: 0 : k as Clearance/Volume (1/h)
ETA3: 0 : Central volume (L)
ETA4: 0 : k21 (1/h)
ETA5: 0 : k12 (1/h)

$PARAM @annotated @covariates
D_WTB : 70 : Body weight (kg) at baseline
EVENING : 1 : evening dose = 1, morning dose = 0
SEX : 0 : 0:Male, 1:Female
AGE : 69 : Age in years at baseline
D_CCRCLB : 70 : Creatinine clearance CCG 

$CMT @annotated
GUT  : Dosing compartment  (mg)[ADM]
CENT  : Central compartment (mg/L)[OBS]
PERIPH: Peripheral compartment V3 ()

$MAIN
double V2 = TVV2*(pow(D_WTB/65,0.806))*(1−0.0217) * exp(ETA3 + ETA(3)) ;
double KA = TVKA * pow((1-0.434), EVENING) * exp(ETA1 + ETA(1)) ; //* pow(0.201, MORNING); // * pow(0.201, MORNING)
double CL = ( (((TVCLR*(D_CCRCLB/70)) + (TVCLNR* (pow(AGE/69, −0.429))*(1+SEX*(−0.215))) )*(1−0.149)) ); // * pow(BW / 70, 1.2) ;
double Q = TVQ  ;
double V3 = TVV3 ;
double k10 = CL/V2 * exp(ETA2 + ETA(2));
double k21 = Q/V3* exp(ETA4 + ETA(4));
double k12 = Q/V2 * exp(ETA5 + ETA(5));

$OMEGA 0.271 0.101 0.0287 0.243 1.49
$SIGMA @labels PROP ADD
0.00 // proportional
0.3 // additive

$ODE
dxdt_GUT = -KA*GUT;
dxdt_CENT = KA*GUT -k12*CENT+k21*PERIPH - k10*CENT;
dxdt_PERIPH = k12*CENT-k21*PERIPH;

$TABLE
double IPRED = CENT/V2;
double DV = (CENT/V2) *(1+PROP)+ADD;

$CAPTURE DV CL KA V2 Q V3
"

grafik

FelicienLL commented 3 years ago

In Cirincione et al, concentrations were log-transformed before analysis, thus the error is not additive, but log-additive or "exponential". You can either tweak your model and your data to work on log-transformed concentrations, or (I would recommend) just define DV and $SIGMA as such in the model:

$SIGMA
0
0.09 // (sigma square = 0.3^2)
$TABLE
double DV = (CENT/V2) * exp(EPS(2));

and mapbayr will automatically log-transform the data in the background for you. Now here is what I get:

Model:  NONMEN_model 
ID : 1  individual(s).
OBS: 2  observation(s).
ETA: 5  parameter(s) to estimate.

Estimates: 
  ID       ETA1       ETA2       ETA3      ETA4       ETA5
1  1 -0.1125337 -0.7029937 0.02217481 0.1390427 -0.4949399

Output (11 lines): 
  ID time evid cmt amt mdv DV      IPRED       PRED       CL        KA       V2    Q   V3 D_WTB EVENING SEX AGE
1  1    0    1   1 2.5   1 NA 0.00000000 0.00000000 2.630706 0.4208702 35.23717 1.62 26.7    80       0   1  76
2  1   12    1   1 2.5   1 NA 0.04192383 0.02835984 2.630706 0.2382125 35.23717 1.62 26.7    80       1   1  76
3  1   24    1   1 2.5   1 NA 0.06913702 0.03845168 2.630706 0.4208702 35.23717 1.62 26.7    80       0   1  76
4  1   36    1   1 2.5   1 NA 0.08680916 0.04695545 2.630706 0.2382125 35.23717 1.62 26.7    80       1   1  76
5  1   48    1   1 2.5   1 NA 0.10108846 0.05007575 2.630706 0.4208702 35.23717 1.62 26.7    80       0   1  76
6  1   60    1   1 2.5   1 NA 0.11047949 0.05484438 2.630706 0.2382125 35.23717 1.62 26.7    80       1   1  76

mapbayr90

andreasmeid commented 3 years ago

Wonderful, Félicien, many thanks! That works perfectly

andreasmeid commented 3 years ago

Close with remark that residual variability needs to be properly defined :-)

FelicienLL / mapbayr

DV missing in output and far too small iiv #91