Download data set only | Datxplore project files | Monolix project files | Mlxplore project files | Simulx scripts
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 1: Introduction
- Part 2: Data visualization with Datxplore
- Part 3: Model development with Monolix
- Part 4: Model exploration with Mlxplore
- Part 5: Simulations for individualized dosing with Simulx and Monolix
Part 3: Model development with Monolix
- Set up of a Monolix project: definition of data, structural and statistical model for the initial model
- Running tasks on the initial model
- Estimation of the population parameter
- Estimation of the standard errors of the parameters, via the Fisher Information Matrix
- Estimation of the individual parameters
- Estimation of the log-likelihood
- Generation of the diagnostic plots
- Model assessment using the diagnostic plots
- Covariate search
- Test of a two compartment model
- Identification of the observation error model
- Convergence assessment
- Overview
Set up of a Monolix project: definition of data, structural and statistical model for the initial model
The setup of a Monolix project requires to define the data set, and the model, both the structural part and the statistical part, including the residual error model, the parameter distributions (variance-covariance matrix) and the covariate model. We will setup a first model and improve it step-wise.
We first create a new project in Monolix and load the data, assigning the columns as we have done in Datxplore. After having clicked “Next”, the user is prompted to choose a structural model.
As a first guess, we try a one-compartment model with bolus administration and linear elimination. This model is available in the Monolix PK library, which is accessed when clicking on “Load from Library > PK”. A (V,Cl) and a (V,k) parameterization are available. We choose the V,k one:
Selecting a model brings us to the next tab “Check initial estimates”, where the user can choose initial parameter values. Although Monolix is pretty robust with respect to initial parameter values, it is good practice to set reasonable initial values as this will speed up the convergence. The simulations with the population parameters are displayed in red on top of the data for each individual in blue. After having increased the number of points per curve in the bottom right corner (“Discretization”) to smooth the displayed curves (important when the Cmax has to be captured accurately), one can play with the values of V and k. In the present case, because the data is sparse, it is not easy to find good initial values. We thus keep the default V_pop=1 and k_pop=1, and click “set as initial values”.
The next tab permits to setup the statistical part of the model, and run estimation tasks.
For the residual error model, we keep the default combined1 model with an additive and a proportional term, and we will see later on if one of the two terms can be dropped.
For the parameter distributions, we keep the default: log-normal distributions for V et k, and no correlations. To start, we also neglect the covariates.
Now that we have specified the data and the model, we save the project as “project_01” (menu: Project > Save). The result text files of all tasks run afterwards will be saved in a folder next to the project file.
All Monolix project files can be downloaded at the end of this webpage.
Running tasks
A single task can be run by clicking on the corresponding button. To run several tasks in a row, tick the box in the upper right corner of the task you want to run and click the green run arrow button. After each task, the estimated values are displayed in the “Results” tab and outputted as text files in the result folder next to the Monolix project.
Estimation of the population parameters
To run the maximum likelihood estimation of the population parameters, click on the “Population parameters” button in the “Tasks” frame. In the pop-up window, you can follow the progression of the parameter estimates over the iterations of the stochastic approximation expectation-maximization (SAEM) algorithm. This algorithm has been shown to be extremely efficient for a wide variety of possibly complex models. In a first phase (before the red line), the algorithm converges to a neighborhood of the solution, while in the second phase (after the red line) it converges to the local minimum. By default, automatic termination criteria are used. In the current case, even of the initial parameter values were pretty far, the algorithm shows a good convergence:
The estimated population parameter values can be accessed by clicking on the “Results” frame.
Estimation of the standard errors, via the Fisher Information Matrix
To calculate standard errors of the population parameter estimates, one can calculate the Fisher Information Matrix (FIM). Two options are proposed: via linearization (toggle “Use linearization method” on) or via stochastic approximation (toggle off). The linearization method is faster but the stochastic approximation is usually more precise. As our model runs fast, we can afford to select the stochastic approximation option and click on the icon to calculate the FIM. After the calculation, the standard errors are added to the population parameter table:
Estimation of the EBEs
The EBEs (empirical bayes estimates) represent the most probable value of the individual parameters (i.e parameter value for each individual). These values are used to calculate the individual predictions that can be compared to the data of each individual to detect possible model mis-specifications of the structural model for instance.
After having run this task, the summary (min, median, max) of the values over all individuals is displayed in the “Results” tab. The estimated value for each individual can be found in the result folder > IndividualParameters > estimatedIndividualParameters.txt file.
Estimation of the conditional distribution
The individual parameters are characterized by a probability distribution that represent the uncertainty of the individual parameter values. The distribution is called conditional distribution and reads 
Estimation of the log-likelihood
For the comparison of this model to more complex models later on, we can estimate the log-likelihood, either via linearization (faster, toggle “use linearization method” on) or via importance sampling (more precise, toggle off). Clicking on the likelihood estimation icon with option “importance sampling” (toggle off), we obtain a -2*log-likelihood of 668, printed in “Results”, together with the AIC and BIC values.
Generation of the diagnostic plots
To assess the quality of the model, we finally can generate built-in graphics. Clicking next to “Plots”, one can see the list of possible graphics, we keep the default list. By clicking on the “Plots” icon, the figures are automatically generated.
Model assessment using the diagnostic plots
In the “individual fits” plot, it seems that the spikes corresponding to the bolus administrations are not well captured, due to insufficient points to draw the prediction. To improve this, go to the wrench button next to Plots, change the “grid” to 1000 and generate the plots again. In the “Individual Fits” plot, you can choose which individual you would like to look at as on the following figure:
One can zoom in to the zone of interest, which is particularly useful for such sparse data sets. On the presented case, we use the X-axis zoom to simplify this action. We here for instance display a zoom on the last time points for the individuals 91 and 92:
The “Observation vs Prediction” plot looks also good:
In the “Distribution of the individual parameters” plot, the theoretical pre-selected parameter distributions (using the estimated population parameters) can be compared to the empirical distribution of individual parameters (histogram). In the “Settings”, if the “conditional distribution” option is chosen (default), the individual parameters for the histogram are drawn from the conditional probability 
|
Conditional distribution
|
Conditional mode
|
|---|
In the “Distribution of the random effects” plot, the standardized random effects plotted as a boxplot are in agreement with the expected standard normal distribution (back horizontal lines represent the median and quartiles). So there is no reason to reject the model for the random effects:
In the “Correlation between random effects” plot, the random effects of the two parameters are plotted against each other. No clear correlation can be seen.
In the “Individual parameters vs covariates” plot, the transformed parameters (log(V) and log(k) as we have chosen log-normal distributions) are plotted against the covariate values. Notice that we use the option “conditional distribution”. This means that instead of using the mode or mean from the post-hoc estimation of each individual parameter, we are sampling from its distribution. This is critical because it is known that the plots using the mode or mean of the post-hoc distribution often show non-existing correlations (see Lavielle, M. & Ribba, B. Pharm Res (2016)).
The covariates figure shows that log(k) decreases with AGE and increases with CLCR. It is possible to display the correlation coefficients in the Settings pop-up window, by selecting “Informations”. Remember that we had noticed with Datxplore that there was a negative correlation between AGE and CLCR. From a biological point of view, a correlation between Tobramycin clearance and creatinine clearance makes sense, as tobramycin is mainly cleared via the kidney.
In summary, the graphics indicate no mis-specification of the structural model or of the residual error model. Yet, the covariates figure suggests to add CLCR as covariate on the elimination k.
Before starting to implement a new model with covariates, we click “Use the last estimates: All” in the Initial estimates tab, to re-use the last estimated parameter values as starting point for SAEM. This button will be disabled (grayed out) as soon as the model is modified.
Covariate search
To add the CLCR covariate on the parameter k, the user can click on the corresponding checkbox in the “Statistical model & Tasks” tab:
The equation corresponding to this click can be displayed by clicking on the “Formula” button: log(k)=log(k_pop) + beta_k_CLCR*CLCR + eta_k. The tick box permits to add the covariate linearly on the transformed parameter. Rewritten, this corresponds to an exponential relationship between the elimination rate k and the covariate CLCR:
$$k_i = k_{pop} e^{\beta_{k,CLCR}\times CLCR_i} e^{\eta_i}$$
with 
In the “Individual parameters vs covariates” plot, the increase of log(k) with CLCR is clearly not linear. So we would like to implement a power law relationship instead:
$$k_i = k_{pop} \left(\frac{CLCR_i}{65}\right)^{\beta_{k,CLCR}} e^{\eta_i}$$
which in the log-transformed space would be written: log(k)=log(k_pop) + beta_k_CLCR*log(CLCR/65) + eta_k
Thus, instead of adding CLCR linearly in the equation, we need to add log(CRCL/65). The next step explains how to create a new covariate logCLCR = log(CLCR/65).
Clicking on the blue button “CONTINUOUS” in the main window, we can define a new covariate called logCLCR equal to log(CLCR/65). The centering permits to ensure that k_pop is representative of an individual with a typical CLCR value (65 in our case). To choose a centering value, hover on the covariate name to see the min, median, mean, and max over the individuals of the data set.
We then add logCLCR as covariate to k, by clicking on the tickbox corresponding to logCLCR and k.
In NONMEM, one would typically write:
TVK = THETA(1)*((CLCR/65)**THETA(2)) K = TVK * EXP(ETA(1))
with THETA(1)=k_pop, THETA(2)=beta and std(ETA(1))= omega_k.
In Monolix, we think in terms of transformed variables:
$$\log(k_i) = \log(k_{pop})+\beta \times \log(CLCR_i/65) + \eta_i$$
but the two models (in Nonmem and Monolix) are the same, corresponding to:
$$k_i = k_{pop} \left(\frac{CLCR}{65}\right)^{\beta_{k,CLCR}} e^{\eta_i}$$
with
We now save the project under project_02, to avoid overwriting the results of the first run. We then select all estimation tasks and click on run to launch the sequence of tasks (Parameter estimation, EBEs, Conditional distribution, Standard errors, Likelihood and Plots).
On the Results frame, we see that the p-value associated with the beta_k_logCLCR parameter is inferior to 1e-10, suggesting that this parameter is significantly different from 0. The effect on CLCR on k is pretty large: for an individual with CLCR=40mL/min (moderate renal impairment) k=0.12 /hr, while for an individual with CLCR=100mL/min (normal renal clearance) k=0.27 /hr. In addition, the standard deviation of k omega_k has decreased from 0.57 to 0.31 as part of the inter-individual variability in k in now explained by the covariate.
In the “Individual parameters vs covariates” plot, we choose to plot the random effects rather than the individual parameters to see if there are remaining correlations (option “Random effects” in “Display”). Note that no correlation is seen anymore between eta_k and AGE, as expected as AGE and CLCR are strongly negatively correlated. It now appears that eta_k is correlated with the weight WT.
|
Individual parameters |
Random effects |
|---|
Depending on the number of individuals, the correlation between random effect and covariate can be more or less significatives. To help the user identifying significative correlations, a Pearson’s correlation test is performed and the p-value is displayed in the Results > Tests tab. Small p-values are highlighted in red/orange/yellow and suggest to test the corresponding parameter-covariate relationship if it also makes biological sense. Note that the correlation tests are performed using the samples from the conditional distribution, which permits to avoid biais due to shrinkage. To learn more about how Monolix circumvents shrinkage, see here.
As the next step, we can add either WT or logWT on log(k). We choose to log-transform and center WT, and add it as a covariate on log(k). We save the project as project_03 and run the sequence of tasks again. beta_k_logWT appears to be significantly different from 0 (p-value below 0.01).
The “Individual parameters vs covariates” plot (but not the statistical tests) then suggests to also add logCLCR on the volume V. We do so and save as project_04. This time the Wald test suggests that logCLCR on log(V) is not significant (p-value=0.29). We thus remove logCLCR on log(V).
Although the “Individual parameters vs covariates” plot doesn’t hint at it, it is common to have the volume V correlated to the weight WT. We try this in project project_05, and the covariate is significant. We thus end with logWT on log(V) and log(k), and logCLCR on log(k). Note that the statistical tests and plots are giving indications about which model to choose, but the final decision lies with the modeler.
The population parameter estimates for this model (project_05) are the following and -2LL is 552.
Test of a two-compartment model
While visualizing the data set with Datxplore, we had noticed that some individuals where hinting at a peripheral compartment. We thus change the structural model to the bolus_2cpt_Vkk12k21 model from the library. As our data is sparse, we choose to remove the variability on k12 and k21, by unselecting the random effects. This mean that all subjects will have the same parameter value for k12 and k21 (no inter-individual variability). We save as project_06.
At this stage, we would like to better understand the impact of extending a (V,k) model to a (V,k,k12,k21) model. We thus open Mlxplore by clicking in the menu on Export > Export to Mlxplore.
Click here to continue to the Mlxplore part of the case study.
Having better understood the influence of the model parameters, we go back to the Monolix interface and set k12 and k21 initial values to the same order of magnitude as k: k12=k21=0.1, and launch the sequence of tasks, including the log-likelihood.
In the Results tab, we see that: (i) the model with two compartments is significantly better (ΔAIC=33 and ΔBIC=28), and (ii) all covariates are still significant.
Could we have seen in the diagnostic plots of the one compartment model that it was mis-specified? Looking at the “Observation. vs Prediction.” figure in log-log scale, we can see that the one-compartment model under-predicts the low concentrations, which is not the case for the two compartment model any more:
|
one-compartment model |
two-compartment model |
|---|
Identification of the observation error model
Until now, we have used a combined1 error model. The constant term “a” is estimated to be 0.00885 in project_06, which is very small compared to the smallest observation (around 0.1). It make thus sense to try a model with a proportional model only (project_07). The error model can be chosen in the “Statistical & Tasks” tab in the frame “Observation model” via the Error model drop down menu. The proportional model appears to be sufficient (ΔAIC=1.47 and ΔBIC=4.04).
Convergence assessment
As a final check, we can do a convergence assessment, to verify that starting from different initial parameters or from different random number seeds leads to the same parameter estimates. For this, we click on “Assessment” (below the run button), select 5 replicates with estimation of the standard errors and log-likelihood. To randomly pick initial values for all parameters, we select “All” under “sensitivity on initial parameters”.
The displayed figures show a good convergence of all replicates to the same minimum. The first plot shows the evolution of the population parameter estimates over the iterations of the SAEM algorithm for each of the 5 replicates. The second plot shows the final estimated parameter value (dot) and the standard errors (horizontal lines) for each parameter and each replicate. The log-likelihood is also displayed.
Final model overview
Our final model is the following:
- 2-compartment model, parameterized as (V, k, k12, k21)
- the residual error model is proportional
- k and V have inter-individual variability. Their distributions are defined by:
$$k_i = k_{pop}(\frac{CLCR}{67.87})^{\beta_{k,CLCR}}(\frac{WT}{66.47})^{\beta_{k,WT}} e^{\eta_{k,i}}~~\textrm{and}~~ V_i = V_{pop}(\frac{WT}{66.47})^{\beta_{V,WT}} e^{\eta_{V,i}}$$
The estimated parameters and their standard errors are:
Click here to continue to Part 5: Simulations for individualized dosing with Simulx and Monolix
