Pharmacometric models simulation using NONMEM, Berkeley Madonna and R

Q-fitter Inc., 6th Floor, 412 Yeoksam-ro, Gangnam-gu, Seoul 06199, Korea.

It is identical to the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/3.0/).

Abstract

In this tutorial, we introduce a differential equation simulation model for use in pharmacometrics involving NONMEM, Berkeley Madonna, and R. We report components of the simulation code and similarities/differences between software, rather than how to use each software. Depending on the purpose of the simulation, an appropriate tool can be selected for effective communication.

Introduction

One of most important steps in pharmacometric (PM) analysis is simulation of various scenarios. Among the many software packages used in PM analysis, NONMEM[1] is still accepted as the gold standard, although the user interface is not as good as other software and it has a steep learning curve.

The purpose of this tutorial is to compare the characteristics of NONMEM, Berkeley Madonna,[2] and R[3] by simulating PM models. It is very important for the modeler to have various software options because the appropriate method can vary depending on who will see the PM analysis results. In this tutorial, we will not introduce the detailed analysis procedure for each program, as our intended audience is those who have a basic understanding of modeling using NONMEM. There are also some excellent tutorials available for PM model simulation using Berkeley Madonna[4] and R.[5,6] The main focus of this manuscript is learning the components of simulation code and the model translation method of this software by working with simple models such as the 2-compartment pharmacokinetic (PK) model and turnover pharmacodynamic (PD) model, as well as dosage regimen adjustment and incorporation of interindividual variability. More complex structural models such as the physiologically-based PK model, target-mediated drug disposition model, complex absorption models, and categorical data simulations will be covered in the next tutorial.

Software

In this tutorial, we will use NONMEM 7.3, Berkeley Madonna 8.3.18, R 3.3.2, and RStudio 1.0.136.[7] NONMEM was developed by ICON Development Solutions (http://www.iconplc.com/technology/products/nonmem/) under a proprietary software license. Berkeley Madonna also requires a license purchase, but most models can be run with the demo version except functions such as ‘save model file’ and ‘export model result.’ R was released under an open-source license.

Example

A simple PK/PD model is used to examine structural model change, dosage regimen adjustment, and incorporation of population variability. The structural model used in the simulation is a 2-compartment disposition PK, inhibition of Kin (synthesis rate of PD) turnover PD model ( Fig. 1 ). Table 1 compares software-specific examples of the simulation code components such as parameters, dosing, sampling, differential equation, integration, and output of the simulation.

Table 1

Example 1: 2-compartment intravenous (IV) bolus PK

NONMEM

In the first example, we assumed a 2-compartment IV bolus PK model result. The parameters were clearance (CL), central volume (Vc), inter-compartmental CL (Q), and peripheral volume (Vp), which were 10 L/hr, 30 L, 10 L/hr, and 50 L, respectively. For the convenience of explanation, there was no inter-individual variability (all OMEGA values were fixed to 0) and residual variability was minimal (additive and proportional error were 0.001 and 0.01). However, because the individual prediction (IPRED), a simulation result of inter-individual variability in the structural model, was used for comparison of the simulation results, residual variability had no effect on simulation performance. In terms of how to perform the simulation, the biggest difference between NONMEM and other software is the necessity of a simulation dataset for simulation in NONMEM. Dosage regimens or sampling points were set in the simulation dataset. The NONMEM code is presented in the following section. Although a 2-compartment IV bolus model can be implemented using ADVAN3 in NONMEM, ADVAN13 (or ADVAN 6) was used because comparison of differential equation models was the main purpose of this tutorial. $MODEL defines the compartments of the model, and PK parameters are defined using the THETA (n) and ETA(n) in $PK. In $DES, differential equations are constructed with appropriate rules. In $ERROR, the gap between individual prediction (IPRED) and observation can be set with residual error. Finally, the variables needed for comparison of the simulation results were stored in a table file for post-processing with R.

$PROB Ex1_2-compartment IV bolus

$DATA simdata_Ex1.csv IGNORE=@

$INPUT ID TIME DV MDV AMT CMT

$SUBROUTINE ADVAN13 TOL=9

DADT(1) = -K10*A(1) - K12*A(1) +K21*A(2)

DADT(2) = K12*A(1) - K21*A(2)

$SIMULATION (20170401) ONLYSIM

$TABLE ID TIME AMT DV MDV IPRED PRED ONEHEADER NOPRINT NOAPPEND

Berkeley Madonna

METHOD, STARTTIME and STOPTIME specify the integration method (RK4; Runge–Kutta 4), start time and stop time of the integration interval. The default integration method is RK4, but other methods such as Euler, RK2, Auto, and Stiff are also available; please refer to the ‘ODE SOLVERS’ section of the reference.[4] DT and DTOUT define the time intervals to be used in the numerical solution of the differential equation system and output time interval, respectively. The default value of DTOUT is zero, and if set to zero, the time interval of DT is applied to the simulation output. Note that the smaller the DT, the longer the time required for the simulation. However, if the DT is too large, the simulation result may change. Therefore, it is necessary to check the difference in simulation results according to the DT. The PK parameter setting is the same as in NONMEM, but direct values are used instead of terms such as THETA and ETA. The differential equations are nearly the same. There are several syntax options for defining differential equations, but the following two methods are the most commonly used.

In NONMEM, since the default value of the compartment at time 0 is 0, the syntax error does not occur even if the compartment value is not specified; however, in Berkeley Madonna, an initial compartment value must be specified even if it is 0. In the case of IV bolus, the initial dose of the central compartment was set as the dose administered. In NONMEM, the scale parameter is reflected and the simulated data represent concentration instead of the amount. However, the Berkeley Madonna does not have such a function, so users must directly define the concentration using the amount and appropriate parameters such as the volume of distribution of central compartment. A ‘DISPLAY’ statement is not mandatory, but you can control which variables are included in the built-in plot window. To specify the variable name to be included in the slider, write DISPLAY again. The sliders window provides a convenient way of changing parameters and automatically running the model. It is common practice to export the data table as a text or comma-separated value (csv) file. In the Berkeley Madonna graph window, clicking on the icon displays two squares partially overlaid. The resulting table gives the data shown in the visual display. The menu “File Save Table As” now allows storage of data as a text or csv file.

R

To solve the differential equation using R, a differential equation solver is required. The default solver provided in R is the ‘deSolve’ package. However, pharmacometric model simulation using the ‘deSolve’ package requires much more complex R coding for dosing regimens or sampling design. Several packages such as ‘PKPDsim,’ ‘RxODE,’ and ‘mrgsolve’ have been developed to overcome these drawbacks. In this tutorial, we used the ‘RxODE’ package because of the ease of comparison with other codes and the author's experience of using this package for simulation in R. A detailed description of ‘RxODE’ and R codes with various examples are provided in the supplementary material in the reference [6]. To use the ‘RxODE’ package, you need to install it by running the following script in the console. Windows users will need to have the appropriate version of Rtools installed. ‘RxODE’ is available on Github at (https://github.com/hallowkm/RxODE) and detailed instructions are provided on this site.

The code structure and the writing of differential equations are very similar to those of the Berkeley Madonna. A compilation manager translates the ODE model into C, compiles it, and dynamically loads the object code into R for improved computational efficiency.[8] Each equation is separated by a semicolon. Dosing and sampling design can be set very intuitively using ‘eventTable.’ An event table object facilitates the specification of complex dosing regimens and sampling schedules. Detailed examples and methods are described in the reference.[6]

Post-processing using R

The output files obtained from NONMEM and Berkeley Madonna were stored as ‘Ex1_NM’ and ‘Ex1_BM.csv,’ respectively, and these files were imported by R for post-processing. This post-processed output was used for plotting. (See supplements).

Example 2: 2-compartment IV infusion PK and adjustment for dosage regimen

NONMEM

IV infusion can be controlled by the infusion rate and duration. To control the infusion rate in NONMEM, users need to enter the amount to be infused per unit time into the RATE column of the dataset. In Example 2, infusion time is automatically defined as the amount / rate (100 mg / 100 mg/hr = 1 hr). Repeated dosing and dosing intervals should be adjusted in the simulation dataset. The dosage regimen can be controlled using ADDL and II columns. In Example 2, we will show an example of dosing every 12 hr for 3 days. To set up this dosage regimen, set ADDL as 5 and II as 12, and sampling time points should be 0–72 hr. There is no change in the control file except that ADDL and II columns are added to $ INPUT.

timepoint 〈- seq(0, 72, 0.1)