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 mrg*solve* 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 $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.

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

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

RyanHi! 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

Sinziana CristeaPost authorHi, 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.

Ryanhi! Hi! Sinziana, Thank you very much for the comment. I greatly appreciate your feedback.