From publication to simulation: extracting information from literature models – Amikacin case study

When you are starting to learn about the basics of modeling & simulation, or when you are planning to start a new study with an already existing compound, it could be useful to reproduce a model you found in literature and play around with the parameters to get model informed predictions.

This post is a case study on how to extract all the information needed to build a population model with the information from a published paper. If something is not clear on the use of mrgsolve or Shiny, have a look at the previous posts on Part 1 of “Creating a simple pharmacometric Shiny application with mrgsolve in R”.

Objectives

  • Identify the structure of the population PK model and the estimated model parameters
  • Identify how we can add covariates to model parameters 
  • Code a model from scratch in mrgsolve
  • Perform PK simulations of the literature model using mrgsolve

Extracting all information from the publication

The publication that we will use for this post is a publication by De Cock et. al. where they built a population PK model to describe the disposition of amikacin (an antibiotic) in preterm neonates after an IV infusion.

Structural model

From the results section of the paper we learn that they identified a two-compartment model to fit their data best. To increase model stability they chose to estimate the inter-compartmental clearance (Q) as a fraction (fr) of clearance (CL):

 Q = fr × CL

Additionally, they set the peripheral volume of distribution (Vp) to be equal to the central volume of distribution (Vd).

Vp = Vd 

A schematic outline of the commonly used abbreviations and the calculation of rate constants for a 2-compartment model can be found in the figure below.

Q = inter-compartmental clearance


Covariates

Pharmacokinetic models developed in the pediatric population commonly use covariates to account for the high levels of inter-individual variability (IIV) existing in the population. A covariate is an explanatory variable of this IIV. E.g. weight may be explaining the variability in the distribution volume.

In this publication, multiple covariate relationships were established:

Clearance

They quantified a significant improvement in the model fit if birth weight (bBW) and postnatal age (PNA) were implemented as covariates on CL, following an exponential and a linear relationship, respectively.

Furthermore, some of the pediatric patients received ibuprofen (a non-steroidal anti-inflammatory drug [NSAID]) as treatment, together with amikacin. This ibuprofen co-administration was found to also impact the CL of amikacin.

Volume(s) of distribution

Current weight (cBW) was implemented as an exponential covariate on the Vd (and on the Vp since the volumes of both compartments were assumed to be equal).

Retrieving the parameter estimates

All the estimated parameters were presented in Table II of the publication, see a copy of the first columns below. We see the fixed effects (used for the structural model and covariate relationships), the interindividual variability (indicating significant variability on the CL parameter), and a proportional and additive residual error component.

The authors provided us with the simple model estimates, without any covariates, and the final PK model estimates. This table clearly shows that the individual CL of a patient is affected by the bBW, PNA and Ibuprofen administration.


In the parameter column, you can observe that the covariate relationships are related to the median covariate values of the preterm neonates population, e.g. bBW/median. The required information on these medians in this population can be found (in this case) in the footnotes of the table.


Building the model in mrgsolve

Now, we have all the information needed to recreate this model! If you are new with using mrgsolve, check out Part 1 of “Creating a simple pharmacometric Shiny application with mrgsolve in R” which will make it much easier to understand the next part.

First, let’s think of which parameters we need to include in the model file (commonly saved with the .cpp extension) to reproduce this published model for amikacin. In order to do this, we need to name all the parameters from the table and initialize them in the model file.

Below, I altered the published Table II with red notes to show where you can find the values of the parameters defining the covariate relationships and matched them to the variable names I will use in the model code.

EBW = Exponent for birth weight
SPNA = Slope PNA
IBU = Impact of ibuprofen co-administration
ECW = Exponent current weight

Now we can put the final parameter estimates in the mrgsolve model file:

$PARAM

// Structural model parameters
TVCL = 0.0493, TVVC = 0.833, TVVP1 = 0.833, TVQFRAC = 0.415

// covariate relationship parameters
EBW = 1.134, SPNA = 0.213, IBU = 0.838, ECW = 0.919

// initial covariate values 
BW = 1750, CW = 1760, PNA = 2, NSAID = 0 

In this case, the initial covariate values were initialized to the median values. These values will be overwritten from the R-script.

In the $CMT section, we say that we have 2 compartments, a central (CENT) and a peripheral (P1) compartment.

$CMT CENT P1 

After which we need to code every covariate relationship in the model. This can be quite tricky and it is good to pay attention to all the different brackets and interactions here:

$MAIN

double CL = (TVCL * pow(BW/1750, EBW) * (1 + SPNA*(PNA/2))) * exp(ETA(1));
double VC = (TVVC * pow(CW/1760,ECW)) * exp(ETA(2));
double VP1 = VC;
double Q1 = (TVQFRAC * TVCL * pow(BW/1750, EBW) * (1 + SPNA*(PNA/2)))*exp(ETA(3)); 

The function pow(number,exponent) is used to code a exponential function in mrgsolve.

The ibuprofen administration is a dichotomous covariate. If the patient receives ibuprofen as co-medication then the NSAID variable takes a value of 1 and CL and Q are multiplied by IBU (because they are both dependent on each other):

if(NSAID == 1) {
CL = CL * IBU;
Q1 = Q1 * IBU;
} 

In the next part we calculate the elimination rate constant (k10) and the distribution rates (k12 & k21), we initialize the $OMEGA and $SIGMA on 0 (a simulation without variability), define our differential equations and define which output we want to capture in a $TABLE:

double k10 = CL/VC;

double k12 = Q1/VC;
double k21 = Q1/VP1;

$OMEGA 0 0 0

$SIGMA @labels PROP ADD
0 0 

$ODE
dxdt_CENT = -k12*CENT+k21*P1 – k10*CENT;
dxdt_P1 = k12*CENT-k21*P1;

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

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

$CAPTURE IPRED DV 

We now save this as new file in the same directory as our R script (amik_popPK.cpp).

Note: When copying text directly from the website, a minus sign can sometimes be converted to a different symbol. Strange error messages in R can therefore be fixed by retypying the minus sign in the .cpp file.

Mrgsolve simulation in R and comparison with literature ( "sanity check")

When you recreate a model from literature, it is of uttermost importance to see whether you can simulate the same profiles as in the article. This will indicate whether your code is correct but also if the authors reported the correct parameters! Sometimes errors can originate from the parameter units (L/h versus L/min), which can quickly result in simulations in the wrong order of magnitude... We will use the bottom row of Figure 4 to compare our model simulations with.

We start with a coding block in which we specify some information in the code, load the required libraries (make sure they are installed), and set the covariates of an individual we would like to simulate.

##################################################################################
############################ CODING BLOCK ########################################
##################################################################################
# Run-Time Environment: R version 3.5.1
# Author: S. Cristea
# Short title: Pharmacokinetic simulations in R for amikacin
# Version: V.1.0
##################################################################################
##################################################################################
###
rm(list = ls(all.names = TRUE))
###

version <- 1.0

###############################################
############ Load libraries
library(mrgsolve)
library(dplyr)
library(ggplot2)

##################################
## Covariate estimates 
bw <- 3520 # Birth weight (g)
pna <- 5 # Postnatal age (days)
cw <- 3800 # Current weight (g)
nsaid <- 0 # Ibuprofen administration (1 - yes, 0 - no) 

In this scenario, we will simulate the profile of an individual patient with a birt weight of 3520, a PNA of 5 days, a current weight of 3800 and no co-administration ibuprofen.

Then, we need to include the dosing information. In pediatric patients doses are commonly given based on their current weight in mg/kg which requires a few extra lines of code compared to a fixed dose. Since the dose will change between individuals, also the rate parameter needs to be different for each patient. However, the duration of the infusion will remain the same and we will therefore define a duration parameter.

###############################################
### Dosing

inf <- T
dose <- 12 * (cw/1000) # mg/kg
ii1 <- 24 # dosing interval in hours

## Infusion duration

d2 <- 20 # min
dur <- d2 / 60 # duration in hours 

Now, it is time to insert the information on our simulation. For how long do we want to simulate? How many individuals do we want to simulate? What percentile intervals do we want to show in our figure?

And we need to create a dosing dataframe. The calculation of the rate makes sure that, regardless of the dose, the duration is equal to the pre-specified variable.

###############################
## Simulation info

nsamples <- 1 ### Number of simulated individuals
sim_time <- 120 ## Time of simulation

# set probabilities of ribbon in figure
minprobs=0.1
maxprobs=0.9

###################################################
############# Set Dosing objects

if(inf){
## IV infusion
Administration <- as.data.frame(ev(ID=1:nsamples,ii=ii1, cmt=1, addl=9999, amt=dose, rate = dose/dur,time=0))
}else{
## IV BOLUS
Administration <- as.data.frame(ev(ID=1:nsamples,ii=ii1, cmt=1, addl=9999, amt=dose, rate = 0,time=0))
}

## Sort by ID and time
data <- Administration[order(Administration$ID,Administration$time),]

We now have all the required objects in R to start our simulation with mrgsolve. We load in the .cpp file, attach the individual covariates to the dataset, and start the simulation. We then immediately calculate the percentiles of the population.

######################################################################
### Load in model for mrgsolve
mod <- mread_cache("amik_popPK")

## Set covariates in dataset
data$BW <- bw
data$PNA <- pna
data$CW <- cw
data$NSAID <- nsaid

#################################################################
###### Perform simulation
out <- mod %>%
data_set(data) %>%
mrgsim(end=sim_time,delta=sim_time/100, obsonly=TRUE) # Simulate 100 observations, independent of the total simulated time

### Output of simulation to dataframe
df <- as.data.frame(out)

#################################################################
## Calculate summary statistics

sum_stat <- df %>%
group_by(time) %>%
summarise(Median_C=median(DV),
Low_percentile=quantile(DV,probs=minprobs),
High_percentile=quantile(DV,probs=maxprobs)
) 

After running this simulation, you get a dataset with your results (sum_stat) which we can use as input for the final figure with the simulated amikacin PK profile and the specified percentiles. We also add the therapeutic window of amikacin to the plot to judge our simulation.

################################### Graphical output


plot_pk <- ggplot(sum_stat, aes(x=time,y=Median_C)) +
  
  ## Add ribbon for variability
  geom_ribbon(aes(ymin=Low_percentile, ymax=High_percentile, x=time), alpha = 0.15, linetype=0)+
  
  ## Add median line
  geom_line(size=2) +
  
  # scale_y_log10()+
  ## Add therapeutic target:
  geom_abline(slope = 0, intercept = 35, col = "red", linetype = "dashed")+
  geom_abline(slope = 0, intercept = 24, col = "green", linetype = "dashed")+
  geom_abline(slope = 0, intercept = 3, col = "red", linetype = "dashed")+
  geom_abline(slope = 0, intercept = 1.5, col = "green", linetype = "dashed")+
  
  # Set axis and theme
  scale_x_continuous("Time (h)",breaks=seq(0,125,24))+
  scale_y_continuous("Concentration (mg/L)",limits=c(0,40))+
  theme_bw()+
  
  # Remove legend
  theme(legend.position="none") 

## Save the plot
png("Simulated_Amikacin.png",width = 10, height = 7.5, units='in',res=600)

plot_pk

dev.off()

Done!

Let’s reproduce the model-informed dosing recommendation from Table IV for patient 3 displayed in Figure 4 – bottom row. This patient should have a bBW = 3520 g. However, since other details on the covariates of this individual needed are missing, we’ll assume cBW = 3800 g and a PNA = 5 days. For this group of cBW > 2800g we are advised to administrate a dose of 12 mg/kg every 24 hours as an IV infusion with the duration of 20-30 min for 125 hours.

Top: Dosing simulation for Patient 3.
Bottom: mrgsolve model simulation.

As you can see, our model prediction closely resembles the information published!

We can now use this model to study different dosing regimens, study the effect of different covariate distributions, assses the variability in the population and even perform clinical trial simulations! Good luck!

GUEST AUTHOR

THIS POST WAS WRITTEN BY Sinziana Cristea FROM Leiden University.
She works on predictive models to study maturation of renal elimination processes, ultimately to guide dosing of renally excreted drugs in children.

COMMENT

Any suggestions or typo’s? Leave a comment or contact me at info@pmxsolutions.com!

2 thoughts on “From publication to simulation: extracting information from literature models – Amikacin case study

  1. Ryan Reply

    Hi! Sinziana, thank you for the tutorial. i have used your approach for a publication. unfortunately i was not able to simulate the numbers in the publication. they are off by ~ 1000 fold. makes me wonder if there is any issue with the units. including details of the publication and my code. any help with the code is greatly appreciated.
    Publication: https://journals.lww.com/drug-monitoring/fulltext/2008/08000/Population_Pharmacokinetics_of_Intravenous,.13.aspx (PMID: 18641540 DOI: 10.1097/FTD.0b013e3181816214)

    ###############################################
    ############ Load libraries
    library(mrgsolve)
    library(dplyr)
    library(ggplot2)

    ##########
    ## Covariate estimates

    wt<- 80

    lbw <- 53

    ###############################################
    ### Dosing

    IM <- F

    dose <- 2 # mg

    #################################################
    ## Insert residual error

    sigmaprop <- 0.10 # Proportional error
    sigmaadd <- 0.001 # Additive error

    ###############################3
    ## Simulation info

    nsamples <- 25 ### Number of simulated individuals
    sim_time <- 10 ## Time of simulation

    # set probabilities of ribbon in figure
    minprobs=0.1
    maxprobs=0.9

    ###################################################
    ############# Set Dosing objects

    if(IM){
    ## IM dose
    Administration <- as.data.frame(ev(ID=1:nsamples, cmt=1, amt=dose, rate = 0,time=0))
    }else{
    ## IV BOLUS
    Administration <- as.data.frame(ev(ID=1:nsamples, cmt=2, amt=dose, rate = 0,time=0))
    }

    ## Sort by ID and time
    data <- Administration[order(Administration$ID,Administration$time),]

    ######################################################################
    ### Load in model for mrgsolve
    mod <- mread("Naloxone")

    ### set covariates in dataset

    data$WT <- wt
    data$LBW <- lbw

    #################################################################
    ###### Perform simulation
    out %
    data_set(data) %>%
    mrgsim(end=sim_time,delta=sim_time/100, obsonly=TRUE)

    ### Output of simulation to dataframe
    df <- as.data.frame(out)

    #################################################################
    ## Calculate summary statistics

    #sum_stat %
    # group_by(time) %>%
    #summarise(Median_C=median(DV),
    # Low_percentile=quantile(DV,probs=minprobs),
    # High_percentile=quantile(DV,probs=maxprobs)
    #)

    sum_stat %
    group_by(time) %>%
    summarise(Median_C=median(DV),
    Low_percentile=quantile(DV,probs=minprobs),
    High_percentile=quantile(DV,probs=maxprobs)
    )

    ################################### Graphical output

    ggplot(sum_stat, aes(x=time,y=Median_C)) +

    ## Add ribbon for variability
    geom_ribbon(aes(ymin=Low_percentile, ymax=High_percentile, x=time), alpha = 0.15, linetype=0)+

    ## Add median line
    geom_line(size=2) +

    # scale_y_log10()+

    # Set axis and theme
    ylab(paste(“Concentration”,sep=””))+
    xlab(“Time after dose (h)”)+
    theme_bw()+

    # Remove legend
    theme(legend.position=”none”)

    #############################

    Drug model file (Naloxone.cpp)

    $PARAM
    //Structural model parameters
    TVKA=0.65, TVCL=91, TVVC=2.87,TVQ1=5.66,TVVP1=1.49,TVQ2=29.8, TVVP2=33.6
    //initial covariate values
    lbw=53, wt=80
    $CMT
    EV CENT P1 P2
    $MAIN
    double KA = TVKA;
    double CL = (TVCL*pow(lbw/70, 0.75))* exp(ETA(1));
    double VC = (TVVC*(wt/70))*exp(ETA(2));
    double Q1 = TVQ1;
    double VP1 = TVVP1;
    double Q2 = TVQ2;
    double VP2= TVVP2;
    double k10=CL/VC;
    double k12=Q1/VC;
    double k13=Q2/VC;
    double k21=Q1/VP1;
    double k31=Q2/VP2;
    $OMEGA 0.00581 0.25
    $SIGMA @labels PROP ADD
    0 0
    $ODE
    dxdt_EV = -KA*EV;
    dxdt_CENT = KA*EV+k21*P1+P2*k31-(k12+k13+k10)*CENT;
    dxdt_P1 = CENT*k12 – P1*k21;
    dxdt_P2 = CENT*k13 – P2*k31;
    $TABLE
    double IPRED = CENT/VC;
    double DV = IPRED*(1+PROP)+ADD;
    while(DV < 0) {
    simeps();
    DV = IPRED*(1+PROP)+ADD;
    }
    $CAPTURE IPRED DV

  2. Sinziana Cristea Post authorReply

    Hi, Ryan.
    Your intuition was right. I also think it’s a unit problem.
    In your code, you give your dose in mg whereas the concentrations in the published plots are in ng/ml. If you give your dose in micrograms instead (1 mg = 1000 micrograms) then the units would match between your simulation and the publication.

    Overall, I think your code looks fine with one small comment: in “Naloxone.cpp” your covariates (lbw, wt) are in lower case:


    //initial covariate values
    lbw=53, wt=80
    […]
    double CL = (TVCL*pow(lbw/70, 0.75))* exp(ETA(1));
    double VC = (TVVC*(wt/70))*exp(ETA(2));

    But in your R script you create a dataset with LBW and WT in upper case:

    ### set covariates in dataset
    data$WT <- wt data$LBW <- lbw " Therefore, I think the changes you make in your dataset in R won't be applied to your model and won't affect your simulations. You should change them to lower or upper case everywhere. Hope this helps.

Leave a Reply

Your email address will not be published. Required fields are marked *