search

metrumresearchgroup / mrgsolve

Simulate from ODE-based population PK/PD and QSP models in R
https://mrgsolve.org
129 stars 36 forks source link

Preventing Negative Values in Simulation Results when using individual parameters as inputs #1174

Closed XuefenYin closed 6 months ago

XuefenYin commented 6 months ago

Hi,

I am reaching out to seek your guidance on an issue I've encountered with my simulation results. I've been using individual parameters as inputs for my simulations, but I've consistently obtained negative values, which theoretically shouldn't occur. Despite my efforts to ensure that the CENT compartment remains non-negative through the simulation code, the issue persists.

Thanks in advance, Xuefen

$PROB

Model: pk2cmt

$PARAM @annotated KA : 0.02 : Absorption rate constant 1 (1/hour) CL : 50 : Clearance (L/hour) V2 : 347 : Central volume (L) Q : 10 : Inter-compartmental clearance (L/hour) V3 : 70 : Peripheral volume of distribution (L) ALAG1: 1.5 : Lag time for KA1 (hour) TR : 0.02: Logit transformation factor for F1

$CMT @annotated GUT : Dosing compartment (mg) CENT : Central compartment (mg) PERIPH : Peripheral compartment (mg) AUC : AUC (mg*hour/L)

$GLOBAL // visit one time

define S2 (V2/1000)

define CP (CENT/S2)

$MAIN // visit every time

double K = CL/V2; double K23 = Q/V2; double K32 = Q/V3; double F1 = 1 - (TR/(1+TR));

ALAG_GUT = ALAG1;

$ODE

double CENT2 = std::max(0.0, CENT);

dxdt_GUT = -KAGUT; dxdt_CENT = KAGUT - K23CENT2 + K32PERIPH - KCENT2; dxdt_PERIPH = K23CENT2 - K32*PERIPH; dxdt_AUC = CP;

$CAPTURE @annotated CP : Plasma concentration (mass/time)

kylebaron commented 6 months ago

Hi @XuefenYin -

Thanks for your question. There are some situations when the amount in CENT could turn negative and I'd be cautious about truncating CENT in $DES.

Can you share some code that executes a simulation where CENT goes to zero? Ideally this would be done via the reprex package https://reprex.tidyverse.org/ . Otherwise, please post a working example in a code block so we can easily run each others code.

For example

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 @annotated
KA : 0.02 : Absorption rate constant 1 (1/hour)
CL : 50 : Clearance (L/hour)
V2 : 347 : Central volume (L)
Q : 10 : Inter-compartmental clearance (L/hour)
V3 : 70 : Peripheral volume of distribution (L)
ALAG1: 1.5 : Lag time for KA1 (hour)
TR : 0.02: Logit transformation factor for F1

$CMT @annotated
GUT : Dosing compartment (mg)
CENT : Central compartment (mg)
PERIPH : Peripheral compartment (mg)
AUC : AUC (mg*hour/L)

$GLOBAL // visit one time

#define S2 (V2/1000)
#define CP (CENT/S2)

$MAIN // visit every time

double K = CL/V2;
double K23 = Q/V2;
double K32 = Q/V3;
double F1 = 1 - (TR/(1+TR));

ALAG_GUT = ALAG1;

$ODE

double CENT2 = std::max(0.0, CENT);

dxdt_GUT = -KA*GUT;
dxdt_CENT = KA*GUT - K23*CENT2 + K32*PERIPH - K*CENT2;
dxdt_PERIPH = K23*CENT2 - K32*PERIPH;
dxdt_AUC = CP;

$CAPTURE @annotated
CP : Plasma concentration (mass/time)
'

mod <- mcode("foo", code)
#> Building foo ...
#> done.

mrgsim(mod, ev(amt = 1)) %>% summary()
#>        ID         time            GUT              CENT        
#>  Min.   :1   Min.   : 0.00   Min.   :0.0000   Min.   :0.00000  
#>  1st Qu.:1   1st Qu.: 5.25   1st Qu.:0.6805   1st Qu.:0.05324  
#>  Median :1   Median :11.50   Median :0.7711   Median :0.08789  
#>  Mean   :1   Mean   :11.54   Mean   :0.7091   Mean   :0.06907  
#>  3rd Qu.:1   3rd Qu.:17.75   3rd Qu.:0.8737   3rd Qu.:0.09420  
#>  Max.   :1   Max.   :24.00   Max.   :0.9900   Max.   :0.09541  
#>      PERIPH              AUC               CP        
#>  Min.   :0.000000   Min.   :0.0000   Min.   :0.0000  
#>  1st Qu.:0.002706   1st Qu.:0.3256   1st Qu.:0.1534  
#>  Median :0.010235   Median :1.6621   Median :0.2533  
#>  Mean   :0.009187   Mean   :1.9459   Mean   :0.1991  
#>  3rd Qu.:0.015320   3rd Qu.:3.3353   3rd Qu.:0.2715  
#>  Max.   :0.017524   Max.   :5.0394   Max.   :0.2750

Created on 2024-03-06 with reprex v2.0.2

This doesn't have negative values in the output. But your simulation would have it.

XuefenYin commented 6 months ago

@kylebaron Thanks so much for the quick reply. I attached the simulation code below, which I get the negative CENT and CP. I also attached the idata.10ID.csv file, which contains the dosing record and individual parameters. I greatly appreciate your help with this matter.

Best regard, Xuefen

input data

idata.10ID.csv

code starts below

library(tidyverse) library(mrgsolve) library(dplyr)

code <- ' $PARAM @annotated KA : 0.02 : Absorption rate constant 1 (1/hour) CL : 50 : Clearance (L/hour) V2 : 347 : Central volume (L) Q : 10 : Inter-compartmental clearance (L/hour) V3 : 70 : Peripheral volume of distribution (L) ALAG1: 1.5 : Lag time for KA1 (hour) TR : 0.02: Logit transformation factor for F1

$CMT @annotated GUT : Dosing compartment (mg) CENT : Central compartment (mg) PERIPH : Peripheral compartment (mg) AUC : AUC (mg*hour/L)

$GLOBAL // visit one time

define S2 (V2/1000)

define CP (CENT/S2)

$MAIN // visit every time

double K = CL/V2; double K23 = Q/V2; double K32 = Q/V3; double F1 = 1 - (TR/(1+TR));

ALAG_GUT = ALAG1;

$ODE

double CENT2 = std::max(0.0, CENT);

dxdt_GUT = -KAGUT; dxdt_CENT = KAGUT - K23CENT2 + K32PERIPH - KCENT2; dxdt_PERIPH = K23CENT2 - K32*PERIPH; dxdt_AUC = CP;

$CAPTURE @annotated CP : Plasma concentration (mass/time) '

mod <- mcode("foo", code)

idata.10ID<- read.csv("idata.10ID.csv")

out <- mrgsim_df(mod, idata.10ID, recover = c("CL"), carry.out = "EVID", recsort = 3, tgrid = seq(0, 25000, by = 1))

summary(out)

Save the negative CP

CP_NEGATIVE <- filter(out, CP < 0)

Check if CP_NEGATIVE dataframe has any rows

if (nrow(CP_NEGATIVE) > 0) { print(paste("There are", nrow(CP_NEGATIVE), "negative numbers in the dataset."))

negative_indices <- which(out$CP < 0)

the first 10 rows of CP_NEGATIVE for verification

print("First 10 rows with negative CP values:") print(head(CP_NEGATIVE, 10)) } else {

print("There are no negative numbers in the dataset.") }

kylebaron commented 6 months ago

Hi @XuefenYin -

Please see the output below. Those negative numbers are happening after essentially all the mass has left the system. The ode solver can only give so much precision in that case and those numbers are just error around zero.

A couple of things you can try

  1. Make atol smaller (like 1e-10 or 1e-20 or something like that; that will give you additional precision around those numbers and they might stay small and positive; this didn't work for me to get rid of all the negative values
  2. You can stop the simulation earlier
  3. You can just set those numbers to zero after the simulation has ended

You can read more about this on this wkii page: https://github.com/metrumresearchgroup/mrgsolve/wiki/Controlling-unphysical-negative-values which I lifted from sundials documentation; it talks about these unphysiological negative values. This other page on setting tolerances would be relevant too: https://github.com/metrumresearchgroup/mrgsolve/wiki/General-advice-on-choice-of-tolerances

Hope this helps! 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 @annotated
KA : 0.02 : Absorption rate constant 1 (1/hour)
CL : 50 : Clearance (L/hour)
V2 : 347 : Central volume (L)
Q : 10 : Inter-compartmental clearance (L/hour)
V3 : 70 : Peripheral volume of distribution (L)
ALAG1: 1.5 : Lag time for KA1 (hour)
TR : 0.02: Logit transformation factor for F1

$CMT @annotated
GUT : Dosing compartment (mg)
CENT : Central compartment (mg)
PERIPH : Peripheral compartment (mg)
AUC : AUC (mg*hour/L)

$GLOBAL // visit one time

#define S2 (V2/1000)
#define CP (CENT/S2)

$MAIN // visit every time

double K = CL/V2;
double K23 = Q/V2;
double K32 = Q/V3;
double F1 = 1 - (TR/(1+TR));

ALAG_GUT = ALAG1;

$ODE

double CENT2 = std::max(0.0, CENT);

dxdt_GUT = -KA*GUT;
dxdt_CENT = KA*GUT - K23*CENT2 + K32*PERIPH - K*CENT2;
dxdt_PERIPH = K23*CENT2 - K32*PERIPH;
dxdt_AUC = CP;

$CAPTURE @annotated
CP : Plasma concentration (mass/time)
'

mod <- mcode("foo", code)
#> Building foo ...
#> done.

idata <- read.csv("/Users/kyleb/helpr/idata.10ID.csv")

out <- mrgsim(mod,
                 idata,
                 recover = c("CL"),
                 carry.out = "EVID",
                 recsort = 3,
                 atol = 1e-12, 
                 rtol = 1e-10,
                 tgrid = seq(0, 25000, by = 1))

summary(out)
#>        ID              TIME            EVID               GUT         
#>  Min.   : 1.000   Min.   :    0   Min.   :0.000000   Min.   :  0.000  
#>  1st Qu.: 3.000   1st Qu.: 6205   1st Qu.:0.000000   1st Qu.:  0.000  
#>  Median : 6.000   Median :12470   Median :0.000000   Median :  0.000  
#>  Mean   : 5.501   Mean   :12474   Mean   :0.002358   Mean   :  4.157  
#>  3rd Qu.: 8.000   3rd Qu.:18735   3rd Qu.:0.000000   3rd Qu.:  0.000  
#>  Max.   :10.000   Max.   :25000   Max.   :1.000000   Max.   :132.639  
#>       CENT             PERIPH            AUC               CP         
#>  Min.   : 0.0000   Min.   :0.0000   Min.   :     0   Min.   :  0.000  
#>  1st Qu.: 0.0000   1st Qu.:0.0000   1st Qu.: 22910   1st Qu.:  0.000  
#>  Median : 0.0000   Median :0.0000   Median : 39829   Median :  0.000  
#>  Mean   : 0.9633   Mean   :0.2149   Mean   : 70284   Mean   :  3.065  
#>  3rd Qu.: 0.0000   3rd Qu.:0.0000   3rd Qu.: 87350   3rd Qu.:  0.000  
#>  Max.   :42.5620   Max.   :8.4453   Max.   :195449   Max.   :122.849  
#>        CL       
#>  Min.   :21.06  
#>  1st Qu.:34.57  
#>  Median :39.11  
#>  Mean   :37.91  
#>  3rd Qu.:44.19  
#>  Max.   :56.12

plot(out)


CP_NEGATIVE <- filter(out, CP < 0)
if (nrow(CP_NEGATIVE) > 0) {
  print(paste("There are", nrow(CP_NEGATIVE), "negative numbers in the dataset."))
  negative_indices <- which(out$CP < 0)
print("First 10 rows with negative CP values:")
print(head(CP_NEGATIVE, 10))
} else {

  print("There are no negative numbers in the dataset.")
}
#> [1] "There are 214825 negative numbers in the dataset."
#> [1] "First 10 rows with negative CP values:"
#> # A tibble: 10 × 9
#>       ID  TIME  EVID      GUT      CENT   PERIPH    AUC        CP    CL
#>    <dbl> <dbl> <dbl>    <dbl>     <dbl>    <dbl>  <dbl>     <dbl> <dbl>
#>  1     1  5727     0 3.47e-18 -3.15e-17 6.18e-18 87350. -9.06e-17  21.1
#>  2     1  5728     0 3.47e-18 -1.93e-16 6.52e-18 87350. -5.56e-16  21.1
#>  3     1  5729     0 3.47e-18 -3.54e-16 6.86e-18 87350. -1.02e-15  21.1
#>  4     1  5730     0 3.47e-18 -5.16e-16 7.20e-18 87350. -1.49e-15  21.1
#>  5     1  5731     0 3.47e-18 -6.77e-16 7.54e-18 87350. -1.95e-15  21.1
#>  6     1  5732     0 3.46e-18 -8.39e-16 7.88e-18 87350. -2.42e-15  21.1
#>  7     1  5733     0 3.46e-18 -1.00e-15 8.22e-18 87350. -2.88e-15  21.1
#>  8     1  5734     0 3.46e-18 -1.16e-15 8.56e-18 87350. -3.35e-15  21.1
#>  9     1  5735     0 3.46e-18 -1.32e-15 8.90e-18 87350. -3.81e-15  21.1
#> 10     1  5736     0 3.45e-18 -1.48e-15 9.24e-18 87350. -4.28e-15  21.1

Created on 2024-03-06 with reprex v2.0.2

XuefenYin commented 6 months ago

@kylebaron Thank you for your clear explanation and suggested solutions. Your expertise has resolved my concerns regarding simulation inaccuracies. Considering the requirement to cover all patient time points, I am opting for Solution 3. Your guidance is greatly appreciated.