Select Page

# Tobramycin case study – Part 5: Simulations for individualized dosing with Simulx and Monolix

This case study presents the modeling of the tobramycin pharmacokinetics, and the determination of a priori dosing regimens in patients with various degrees of renal function impairment. It takes advantage of the integrated use of Datxplore for data visualization, Mlxplore for model exploration, Monolix for parameter estimation and Simulx for simulations and best dosing regimen determination.

The case study is presented in 5 sequential parts, that we recommend to read in order:

### Part 5: Simulations for individualized dosing with Simulx and Monolix

1. A glance at Simulx possibilities
2. Toxicity and target attainment of the usual treatment
3. A priori determination of the dosing regimen for an individual
• Using the individuals covariates
• Using the individuals parameters obtained from early drug monitoring
4. A note on individual predictions for individualized dosing

### A glance at Simulx possibilities

Now that we have a model, we can use it to do simulations. The saved Monolix project can be reused by the R package Simulx for straightforward simulations. All project components (structural model, number of subjects, covariates, treatments, estimated parameter values and outputs) can either be reused from the project or data set, or be redefined by the user.

After having open R, loaded the mlxR library (containing simulx) via library(mlxR) and set the working space, we are ready to start. For instance, to simulate the project without any change, just type:

res1  <- simulx(project = '../monolixProjects/project_07.mlxtran')

The results can be plotted with:

print(ggplotmlx(data=res1$y) + geom_point(aes(x=time, y=y, colour=id)) + geom_line(aes(x=time, y=y, colour=id)) + scale_x_continuous("Time") + scale_y_continuous("Observed concentration") + theme(legend.position="none")) We obtain the following spaghetti plot: As a second example, we can change the administration schedule to two doses of 80 mg at 0, 8 and 16 hr for all individuals, simulate only 10 individuals (resampled from the data set with their covariates), change the V_pop parameter to a higher value and set denser measurement times. Only a few lines are necessary to define the simulations to run: #====== define project file project.file <- '../monolixProjects/project_07.mlxtran' #====== define administration schedule adm <- list(time=c(0,8,16), amount=80) #====== define number of individuals N <- 10 #====== define observation times and simulation outputs outy <- list(name = 'y', time = c(0,1,2,5,8,9,10,13,16,17,18,21)) outCc <- list(name = 'Cc', time = seq(0,25, by=0.1)) #===== change V_pop parameter value (others are taken from the estimates.txt file) param <- c(V_pop=30) #===== simulate res2 <- simulx(project = projectfile, output = list(outy, outCc), treatment = adm, group = list(size=N), parameter = param) The predicted exact concentrations and observed concentrations are the following: ### Toxicity and target attainment of the usual treatment In the present case, we are interested in optimizing the dosing regimen of Tobramycin. Usually Tobramycin is administrated as an intravenous bolus at 1mg/kg every 8h. Higher clinical efficacy has been associated with peak drug concentration of 6-7 mg/L. Thus our target attainment is set at 6-7 mg/L, while the trough concentration must be below 2mg/L to avoid toxicity. What percentage of the population meets the efficacy target and the lower limit, if Tobramycin is administrated according to the usual dosing regimen? To answer this question, we will simulate a typical population. In a typical adult population, the mean weight of men is 72 kg with a standard deviation of 10 kg, and the mean weight of women is 57±8.6 kg (FDA, “Physiological parameter values for PBPK models”, 1994). The mean creatinine clearance is 101±30.5 mL/min (Kesteloot et al., J Hum Hypertens., Vol.10, No. 4, 1996). The population of our data set is relatively limited and may not represent well the entire population, especially regarding the creatinine clearance. We will thus try two methods: reuse the covariates from the dataset or draw the covariates from the known population distribution. In the Monolix project, doses are given in mg. We now would like to work with doses in mg/kg, which makes the code slightly more complex than the two examples above. To simulate the steady-state peak and trough concentration for 1000 individuals, resampled from the data set, receiving 1 mg/kg every 8 hours: project.file <- '../monolixProjects/project_07.mlxtran' N <- 1000 # output definition: look at 2 cycles after 35 cycles (to be in steady-state) outCc <- list(name = 'Cc', time = seq(279, 295.9, by=0.1)) # defining the repeated administration. # The amount is calculated based on the weight d <- readDatamlx(project = project.file) covariates1 <- d$covariate
dosePerKg <- 1
times <- seq(0,300,by=8)
nbTimes <- length(times)
adm1 <- data.frame(id = rep(1:d$N,each=nbTimes), time = rep(times,d$N),
amount = rep(dosePerKg*covariates1$WT,each=nbTimes)) # simulation of N=1000 individuals res1 <- simulx(project = project.file, output = outCc, treatment = adm1, parameter = d$covariate,
group = list(size=N))

We can then plot the concentration prediction interval with:

pPrct1 <- prctilemlx(res1$Cc, band=list(number=9,level=90)) + scale_x_continuous("Time (hr)") + scale_y_continuous("Concentration (mg/L)", limits = c(0, 9.5)) + ggtitle("with covariates sampled from the data set") print(pPrct1)  We can now do the same, but instead of resampling the covariates from the data set, we draw the covariates from the known distribution of weight and creatinine clearance. We assume that the population is composed of 50% men and 50% women. The code for the covariates definition reads (the rest of the code is similar to the one above):  # covariates from representative distributions covariates2 <- data.frame(id = (1:N), WT = c(rnorm(n=N/2, mean=72, sd=10),rnorm(n=N/2, mean=57, sd=8.6)),#half men, half women CLCR= rnorm(n=N, mean=101, sd=30), AGE = NaN, #unused but needed for consistency SEX = NaN ) #unused but needed for consistency # defining the repeated administration. The amount is calculated based on the weight adm2 <- data.frame(id = rep(1:N,each=nbTimes), time = rep(times,N), amount = rep(dosePerKg*covariates2$WT,each=nbTimes))

# simulate the 1000 individuals
res2 <- simulx(project = project.file,
output = outCc,
parameter = covariates2)

We obtain the following prediction distribution:

The prediction interval obtained from the distribution representative of the whole population is more narrow than that obtained by resampling the data set. This was expected as the data set includes a large proportion of individuals with impaired renal function (low creatinine clearance).

Using the simulations, one can easily calculate that for 24% of the population the default dosing will be toxic, while for 82% of the population it will be ineffective. This calls for an individualization of the dosing regimen, using the individuals covariates.

### A priori determination of the dosing regimen for an individual

#### Using the individuals covariates

We this would like to a priori determine an appropriate dosing regimen for an individual, knowing its covariates values. For instance, we will consider an individual weighting 78 kg, with a creatinine clearance of 30 mL/min, which corresponds to a severe renal impairment. We first consider a typical individual with this weight and creatinine clearance, i.e we neglect the random effects.

How would the typical treatment perform on this individual? We simulate the patient with the following code:

# define project file
project.file <- '../monolixProjects/project_2cpt_final.mlxtran'

# an individual with WT=78 and CLCR=30
weight <- 78
clcr <- 30

# define output
outCc <- list(name = 'Cc', time = seq(279, 295.9, by=0.01))

# define covariates
covariates <- data.frame(id = 1,
WT = weight,
CLCR= clcr,
AGE = NaN, #unused but needed for consistency
SEX = NaN) #unused but needed for consistency

#disable random effects to simulate typical individual with this WT and CLCR
param=c(omega_k=0,omega_V=0)

# define usual treatment
dosePerKg <- 1
interdoseinterval <- 8

# run simulation of one individual
res <- simulx(project = project.file,
output = outCc,
parameter = list(covariates,param))

After plotting, we obtain the following figures that shows that for this individual the treatment is toxic.

We thus would like to adapt the dosing regimen to this individual, given its covariates values. We start by implementing a function that returns the peak and trough concentrations at steady-state, for a given dose per kg and a given inter-dose interval:

getMinMax <- function(projectfile,dosePerKg,interdoseinterval, weight, clcr){

outCc <- list(name = 'Cc', time = seq(35*interdoseinterval-0.05, 35*interdoseinterval+28, by=0.05))
# define covariates
covariates <- data.frame(id = 1,
WT = weight,
CLCR= clcr,
AGE = NaN, #unused but needed for consistency
SEX = NaN) #unused but needed for consistency

# define treatment

#disable random effects
param=c(omega_k=0, omega_V=0)

# run simulation of one individual
N <- 1
res <- simulx(project = projectfile,
output = outCc,
parameter = list(covariates,param))

# set starting time to zero for easier comparison
Conc <- res$Cc Conc$time <- Conc$time - Conc$time[1]

maxi <- Conc[which(diff(sign(diff(Conc$Cc)))==-2)+1,] mini <- Conc[which(diff(sign(diff(Conc$Cc)))==2)+1,]

result <- list(min=mini$Cc[1], max=maxi$Cc[1])
return(result)
}

We then implement a simple optimization algorithm: if the peak is too high or too low, the dose is decreased or increased respectively. If the trough is too high, the inter-dose interval is increased. The code reads:

#initial dosing regimen
dosePerKg <- 1
interdoseinterval <- 8
weight <- 78
clcr <- 30
result <- getMinMax(dosePerKg, interdoseinterval,weight,clcr)

while(result$max < 6 | result$max > 7 | result$min> 2){ if (result$max < 6){
dosePerKg = dosePerKg + 0.1
}else if(result$max > 7){ dosePerKg = dosePerKg - 0.1 }else if(result$min > 2){
interdoseinterval = interdoseinterval + 1
}else{
paste0("Error")
}
result <- getMinMax(dosePerKg, interdoseinterval,weight,clcr)
}

In the video below, one can visualize how the dose and the inter-dose interval evolve from the initial default treatment to a treatment achieving both efficacy and non-toxicity. For this individual, we propose 1.2 mg/kg administrated every 14 hr.

For the optimization above, we have considered a typical individual among the population of people with WT=78 and CLCR=30. We will now verify the chances of the optimized treatment to be safe and efficient, taking into account the random effects on the parameters for this individual.

The code is as follows. The only difference is that we do not disable the random effects anymore and simulate the individual 1000 times to get a prediction distribution.

#output
outCc <- list(name = 'Cc', time = seq(35*interdoseinterval-0.05, 35*interdoseinterval+29.9, by=0.05))

# to get
N <- 100

# define covariates
covariates <- data.frame(id = (1:N),
WT = weight,
CLCR= clcr,
AGE = NaN, #unused but needed for consistency
SEX = NaN) #unused but needed for consistency

# define usual treatment

# run simulation of one individual, parameter k is drawn
#from its distribution, taking into account the covariates
res <- simulx(project = project.file,
output = outCc,
parameter = covariates)

With the optimized treatment for this individual, there are 10% chance that the treatment is toxic and 31% chances that it is not effective. This is an improvement compared to the default treatment (for this individual, 92% chances toxic and 13% chances not effective). Yet if one would have a better estimate of the parameters for this individual instead of only its covariates, one would be able to even better calibrate the dosing schedule. This is possible using therapeutic drug monitoring.

#### Using the individuals parameters, obtained from early drug monitoring

After the first drug injection, it may be useful to monitor the drug concentration in order to adapt the subsequent doses. The goal is to use the measured concentrations to determine the parameters of this individual and use these estimated parameters to simulate the expected concentration range after several doses and adapt the dosing accordingly.

Let’s assume that for our individual (weight=78 kg, creatinine clearance=30 mL/min), we have started with the a priori proposed dose 1.2 mg/kg every 14 hr, and have measured the concentration at a few time points after the first dose:

ID   TIME  DOSE     Y  WT  CLCR
1     0    93.6     .  78   30
1     1    .     3.59  78   30
1     3    .     2.75  78   30
1     5    .     2.07  78   30
1     8    .     1.11  78   30

We can use Monolix to estimate the individual parameters of this individual, fixing the population parameters to those estimated previously:

After having run the estimation of the population parameters (which just prints the parameters to a file, as all population parameters are fixed), we can estimate the mean, mode and std of the conditional distribution  $p(\psi_i|y_i;\hat{\theta})$ (individual button). The results are printed in the file “indiv_parameters.txt” in the result folder. We obtained the following values: k_mean=0.112, k_std=0.024, V_mean=23.8 and V_std=2.26.

The “IndividualFits” graphic shows the prediction using the estimated conditional mean and the drug monitoring data:

We now would like to use the individual parameters to adapt the treatment. In Simulx, it is not possible to directly set the value of the parameters at the individual level when using a project file as input. We will instead work with the structural model file. In a very similar manner as previously, we define a getMinMax function (below) and an optimization loop.  V, and k are to the mean of the conditional distribution calculated with Monolix.

getMinMax_withModel <- function(modelfile, dosePerKg,interdoseinterval, reference, weight,clcr){

outCc <- list(name = 'Cc', time = seq(35*interdoseinterval-0.05, 35*interdoseinterval+28, by=0.05))

# parameter values fixed
parameters <- data.frame(id = 1,
V = 23.8,
k = 0.112,
k12 = 0.0303,
k21 = 0.1)

# define treatment

# run simulation of this individual
res <- simulx(model = modelfile,
output = outCc,
parameter = parameters)

Conc <- res$Cc Conc$time <- Conc$time - Conc$time[1]

maxi <- Conc[which(diff(sign(diff(Conc$Cc)))==-2)+1,] mini <- Conc[which(diff(sign(diff(Conc$Cc)))==2)+1,]

result <- list(min=mini$Cc[1], max=maxi$Cc[1], ref=reference)
return(result)
}

The optimization leads to the following schedule: 1.6 mg/kg every 18hr.

Due to the uncertainty in the individual parameters for this individual, there is also an uncertainty on the simulated concentration.

To verify the chances of the optimized treatment to be safe and efficient, we would like to take into account the distribution of the individual parameter instead of only its mean. The distribution of the individual parameters is a joint distribution, that can for instance be estimated via a Monte-Carlo procedure. As this is not trivial, we will instead assume that the conditional distribution of the individual parameters Vi and ki is a Gaussian, and reuse the standard deviation calculate by Monolix. The code is as follows:

N <- 1000

# parameter values from distribution
parameters <- data.frame(id = (1:N),
V = rnorm(n=N, mean=18.86, sd=1.66),
k = rnorm(n=N, mean=0.156, sd=0.025),
k12 = 0.149,
k21 = 0.275)

# run simulation of N individual, parameter k is drawn from its distribution
outparam <- c("V","k","k12","k21")
res2 <- simulx(model = structuralModelFile,
output = list(outCc,outparam),
parameter = parameters)

# calculate the percentiles and plot
pPrct2 <- prctilemlx(res2\$Cc, band=list(number=9,level=90)) +
scale_x_continuous("Time (hr)") +
scale_y_continuous("Concentration (mg/L)") + ylim(c(0,7.5))

print(pPrct2)

Using the individual parameters to optimize the treatment, there are 8% chances that the treatment is toxic and 37% chances that it is not efficient. In this case, the information provided by the drug monitoring does not permit to improve the treatment design significantly. Another patient for which the drug monitoring brings a lot is presented in the next section.

### A note on individual predictions for individualized dosing

In the previous sections, we have seen three ways of predicting the concentration time course of an individual:

1. If no covariate information is available for this individual, we only know that the concentration time course for the individual will be within the prediction interval of the entire population, which is pretty large (tan color in the figure below).
2. If the covariates of this individual are available, we can use them to get a reduced distribution of the parameters k and V for this individual. The associated concentration prediction (bleu color in the figure) is narrower compared to the first case.
3. If in addition to the covariates, a few measures can be made (drug monitoring), the conditional distribution of parameters k and V (given theses measures and the previously estimated population parameters) can be obtained. This leads to an even more narrow concentration prediction interval (brown in the figure).

Note that for this figure we have used a more typical individual without renal impairment (CLCR=80 mL/min).

This example shows how a  personalized treatment can be realized using population PK modeling and simulation, and how the MonolixSuite software components facilitate an efficient implementation. The modeling/simulation approach permits to precisely assess the trade-off between the prediction precision and the costs of information acquisition. Covariates are usually relatively cheap to measure but in the present case lead to only a moderate reduction of the prediction interval. On the other hand, therapeutic drug monitoring is more expensive but permits a strong narrowing of the concentration prediction interval, which in turn permits a better determination of an effective and non-toxic dosing regimen.