# 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.

#### 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.

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

0 0

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

\$TABLE
double IPRED = CENT/VC;

while(DV < 0) {
simeps();
}

\$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

###############################################
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
}else{
## IV BOLUS
}

## Sort by ID and 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

## 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)) +

geom_ribbon(aes(ymin=Low_percentile, ymax=High_percentile, x=time), alpha = 0.15, linetype=0)+

geom_line(size=2) +

# scale_y_log10()+
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.

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!

COMMENT

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

1. 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)

###############################################
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

###############################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

######################################################################
### Load in model for mrgsolve

### 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)) +

geom_ribbon(aes(ymin=Low_percentile, ymax=High_percentile, x=time), alpha = 0.15, linetype=0)+

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
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;
while(DV < 0) {
simeps();
}
\$CAPTURE IPRED DV

2. 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.