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.

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.

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]

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

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.

In addition to implementing IV infusion, we will look at codes that allow flexible setting of the number of doses and dosing interval in Berkeley Madonna. First, dose, duration of infusion, infusion rate, dosing interval, and number of doses are set to dose, inf_duration, inf_rate, inf_interval, and ndoses. The ‘is_inftime’ variable is an indicator that allocates 1 during the infusion and 0 for after the infusion ends. The ‘dosingperiod’ is also an indicator variable determined by the number of doses (ndoses) and dosing interval (inf_interval): 1 for the period where dosing is going on and 0 for the period without dosing. The actual infusion consists of a combination of three variables, inf_rate, is_inftime, and dosingperiod, which are assigned to the ‘input’ variable. The input is incorporated in the dosing compartment of the differential equation (in this case, A1).

In the R code of Example 2, you only need to change the add. dosing and add.sampling parts of the eventTable as follows:

Implementation of oral absorption is relatively simple. An absorption compartment should be added to $MODEL, a differential equation for first-order absorption in $DES, and a first-order absorption rate constant (KA) parameter are also added to $PK.

In this case, dosing is entered in the absorption compartment (CMT=1) and the observation is made in the central compartment (CMT=2), so the number of CMT columns in the simulation dataset must be adjusted accordingly.

The Berkeley Madonna also requires an additional compartment, differential equation, and KA parameter. Since it is not an infusion, only the indicator of the dosing period is set up for repeated administration. After defining the input variable using the pulse function, we put the input into the A1 compartment. The ‘pulse’ function, which is a Berkeley Madonna built-in function that allows you to set the dose, lag time, and dosing interval, is often used in dosing settings.

The same principle applies to R, and the part of the code that is changed is as follows.

First, the sampling points of the observation compartment in the simulation dataset should be set to PK and PD, respectively, and distinguished in the CMT column. Therefore, PD sampling is added to the existing dataset and the compartment number is set to 4. The R code for dataset generation is as follows.

In the control file, it is necessary to set the initial value of the PD compartment as baseline (A_0(4) =

BASE) in addition to adding the PD compartment to the $MODEL defining the PD parameter and correcting the differential equation. Since PK and PD observations must be computed separately, we set the error models for CMTs of 2 and 4 at $ERROR, respectively.

PD compartment, differential equation, PD parameters, and initial value of the PD compartment are set. The changes compared to Example 3 are as follows.

The same principle applies to R, and the part of the code that is changed is as follows.

Example 5 simulates 1000 individuals with 20% variability in CL, Vc, and Ka with the 2-compartment oral PK model. The only changes are the number of ID (nID) value from 1 to 1000 in the dataset simulation code, and $OMEGA in the NONMEM control file given as 0.04 for CL, Vc, and Ka instead of ‘0 FIX.’

In Berkeley Madonna, array equations are useful tools that can be used for many purposes. One-dimensional arrays (vectors) are defined as follows:

Arrays of differential equations are defined by placing the brackets immediately following the variable name. For example, the following equations define a set of N differential equations:

In addition to the number of ID (N) and standard deviation of the parameter (CLsd, V2sd, KAsd), most of the code has been changed to the array type compared with the Ex3 code.

In R, we need to generate an individual ETA matrix to reflect interindividual variability. In this tutorial, the ‘mvrnorm’ function of the ‘MASS’ package was utilized and an individual parameter matrix was generated. Each individual was simulated by looping through the parameter matrix, using each row as an input for the simulation, and collecting the output of each simulation in an output matrix.