# 2  Model specification

This chapter details the mrgsolve model specification format.

## 2.1 How / where to write a model

There are two ways to write your model:

1. Code in a separate file and source into your R script
2. Code inline as a character string already in your R script

We recommend method 1 (separate file) for any non-trivial modeling work. Method 2 is handy for quickly coding a model and you’ll also see us using that approach frequently when demonstrating how to use mrgsolve.

### 2.1.1 Separate file

For most applications, you will want to put your model code in a standalone file. This file can have any extension, but there is some special behavior when you use the .cpp extension. We’ll show you both here.

#### 2.1.1.1 Using .cpp file extension

Open a text editor and type the model code into a file with name that has the format <model-name>.cpp. This filename format identifies a “name” for your model (<model-name>, the “stem” of the file name). Note: this whole file will be read and parsed, so everything in it must be valid mrgsolve model specification elements.

Use the mread() function to read and parse this file. For the model called mymodel saved in mymodel.cpp (in the current working directory), issue the command:

mod <- mread("mymodel")

mread() returns a model object from which you can simulate.

#### 2.1.1.2 Using any file extension

You can use any file extension to save your mrgsolve model. Rather than using .cpp, you might use .mod. When you are using an extension other than .cpp, pass in the entire file name, both stem and extension. If you have the model code saved in mymodel.mod then you’d call

mod <- mread("mymodel.mod")

#### 2.1.1.3 The model project directory

The second argument to mread() is project. This is the directory where mrgsolve will look for your model file (passed as the first argument to mread()).

If your myproject.mod file is located in the models directory, then you could load that file with

mod <- mread("mymodel.mod", project = "models")

Alternatively, you can just past the entire path to the model file

mod <- mread("models/mymodel.mod")

### 2.1.2 Inline / code

Often it is more convenient to write a model right in your R script. The model might look something like this:

code <- '
$PARAM CL = 1, VC = 20$PKMODEL ncmt=1
'

Here, we created a character vector of length 1 and saved it to the R object called code. The name of this object is irrelevant. But code will be passed into mrgsolve as the model definition. When mrgsolve gets a model like this along with a “name” for the model, mrgsolve will write the code to a file called <model-name>.cpp and read it right back in as if you had typed the code into this file (Section 2.1.1).

To parse and load this model, use the mcode() command:

mod <- mcode("mymodel", code)

mcode() is a convenience wrapper for mread(). mcode writes the code to mymodel.cpp in tempdir(), reads it back in, compiles and loads.

The mcode call is equivalent to:

mod <- mread("mymodel", tempdir(), code)

For help, see ?mread , ?mcode in the R help system after loading mrgsolve.

You can comment out code in the model specification file with two front slashes //

$MAIN double a = 1.234; // double b = 5.678;  This comment sequence is what is used in C++ code, but you can use it anywhere in the mrgsolve model specification file. ## 2.2 Code blocks ### 2.2.1 About code blocks Block identifier Different types of code are organized in the model specification file and separated by block identifiers. There are two ways to formulate block identifiers that can be used in mrgsolve. In the first type, a dollar-sign is placed at the start of the block name $BLOCKNAME
<block-code>

For example, a block of parameters would be

$PARAM CL = 1 The second way to write this is with brackets [ BLOCKNAME ] <block-code> There is no functional difference between the dollar-sign notation and the brackets. When model specification code is saved into a file with a .cpp extension, the code editor may make certain assumptions about formatting or styling the code. Using brackets will most-likely work better with the editor in that case. Block identifiers are case-insensitive so all of these also work $param
CL = 1
[ param ]
CL = 1

Users are free to include block code on the same line as the block identifier, but must include a space after the identifier. For example, the parser will recognize $PARAM CL = 1 but not $PARAMCL=1 as parameters.

Block syntax Different blocks may require different syntax. For example, code written in $PARAM will be parsed by the R parser and will generally need to adhere to R syntax requirements. On the other hand, code in $MAIN, $ODE, and $TABLE will be used to make functions in C++ and therefore will need to be valid C++ code, including terminal ; on each line.

Block options Options may be specified on some code blocks that signal how the code is to be parsed or used in the simulation.

Block options can be boolean in which case you indicate the option name with nothing following. For example, to indicate block = TRUE for an OMEGA matrix, write

$OMEGA @block 1 2 3 Here @block indicates block = TRUE. Starting with mrgsolve 1.0.0, block options can be negated by writing ! between @ and the option name. To request a non-block (diagonal) OMEGA matrix $OMEGA @!block
1 2 3

But note well that diagonal is the default configuration for OMEGA, so it is never necessary to write this.

Other non-boolean block options state the name and then a value to be assigned to that name. For example, if we wanted to assign labels to the ETAs from an OMEGA block, we’d write

$OMEGA @labels ECL EVC 1 2 This essentially sets labels equal to the vector c("ECL", "EVC"). We note two specific boolean options (@object and @as_object) that have similar function across multiple blocks. ### 2.2.2 Programmatic or bulk initialization The following describes syntax for initializing blocks as R objects using R code. This can be helpful, say, when you need to initialize a 50x50 OMEGA matrix or a series of systematically named parameters. We will describe the @object and @as_object options that are available on select blocks. The @object option lets you name an object defined in $ENV to use to instantiate the block data. For example, to specify a series of parameters using the @object option, you’d write

$ENV params <- list(CL = 1, V = 20, KA = 1.2)$PARAM @object params

This tells mrgsolve that there is an object called params that it is in the $ENV environment and mrgsolve will use that to define the names and values. This is a trivial example to show how a simple series of parameters could be defined. However, the intended use for this functionality is to allow efficient creation of a large, systematically-named series of parameters in a large model. The @as_object option is a boolean option that tells mrgsolve that the block code will actually return the object (rather than asking mrgsolve to look in $ENV for the object). The equivalent specification for the block above would be:

$PARAM @as_object list(CL = 1, V = 20, KA = 1.2) The following blocks contain both the @object and @as_object options: • $PARAM
• $THETA • $CMT
• $INT • $OMEGA
• $SIGMA Please see the specific block documentation for more details on the specific type of object that should be returned when this syntax is invoked. ### 2.2.3$PROB

Syntax: text
Multiple allowed: yes
Options: @annotated, @covariates

Use this block to make notes about the model. There are no restrictions on the text that gets entered here. mrgsolve does not routinely process the text in any way, except when rendering the model as a document. Frequently, we write the text here in markdown format so that it will render nicely in the model document. But this is completely optional.

See the annotated model in Section 11.1 for an example.

### 2.2.4$PARAM Syntax: R Multiple allowed: yes Options: @annotated, @covariates, @object, @as_object Define the parameter list in the current model. Parameters are names associated with values that can be used throughout the model. A value must be given for every parameter name. Names (and numbers) of parameters must be set at the time the model is compiled, but parameter values may be updated without re-compiling the model. The @covariates option allows you to tag a collection of parameters as “covariates”. It does not change the functionality of the model or simulation workflow in any way, but allows you to get that list of covariate names out of the model object. Example: [ PARAM ] CL = 1, VC = 20, KA = 1.2 KM = 25, VMAX = 400, FLAG = 1, WT = 80 SEX = 0, N = sqrt(25) Note that the “values” in the parameter list will get evaluated by the R interpreter and the evaluation will happen inside the environment which is defined by the $ENV block. For example

$ENV NWT = 1.2$PARAM
WT = 80 * 2.2
MASS = 600 * MWT

Annotated example:

[ PARAM ] @annotated
CL :   1 : Clearance (L/hr)
VC :  20 : Volume of distribution (L)
KA:  1.2 : Absorption rate constant (1/hr)

See the popex model for an example of using @covariates

.
. Model file: popex.cpp
.
. $PARAM . TVKA = 0.5, TVCL = 1, TVV = 24 . .$PARAM
. @covariates
. WT = 70

then

mod <- modlib("popex")
as.list(mod)$covariates . [1] "WT" Notes: • Multiple blocks are allowed • Values are evaluated by the R interpreter The @object and @as_object options should name or return a named list of parameters. See also: Section 2.2.22 and Section 2.2.5. See ?param in the R help system after loading mrgsolve. ### 2.2.5$FIXED

Syntax: R
Multiple allowed: yes
Options: @annotated, @object, @as_object

Like $PARAM, $FIXED is used to specify name=value pairs. Unlike $PARAM, however, the values associated with names in $FIXED are not able to be updated.

By default, names in $FIXED are associated with their value through a C++ preprocessor #define statement. Usually, $FIXED is only used when there are a very large number of parameters ($$>$$ 100 or 200). When some of these parameters never need to be updated, you can move them to a $FIXED block to get a modest gain in efficiency of the simulation. Items in $FIXED will not be shown when parameters are queried.

Example:

[ PARAM ] CL = 2, VC = 20

[ FIXED ]
g = 9.8

Annotated example:

$FIXED @annotated g : 9.8 : Acceleration due to gravity (m/s^2) See also: Section 2.2.4 and Section 2.2.22. Notes: • Multiple blocks are allowed • Values are evaluated by the R interpreter ### 2.2.6$CMT and $INIT Syntax: text Multiple allowed: yes Options: @annotated, @object, @as_object Declare the names of all compartments in the model. • For $CMT give the names of compartments; initial values are assumed to be 0
• For $INIT give the name and initial value for all compartments Note that both $CMT and $INIT declare compartments, so any compartment name should get declared in either $CMT or $INIT, but never both. Examples: [ CMT ] GUT CENT RESPONSE [ INIT ] GUT = 0, CENT = 0, RESPONSE = 25 Annotated examples: [ CMT ] @annotated GUT : Dosing compartment (mg) CENT : Central PK compartment (mg) RESPONSE : Response $INIT @annotated
GUT      :   0 : Dosing compartment (mg)
CENT     :   0 : Central PK compartment (mg)
RESPONSE :  25 : Response

The @object and @as_object options should name or return a named list of compartments and initial values when used in $INIT and a character vector of compartment names when used in $CMT.

See ?init in the R help system after loading mrgsolve.

### 2.2.7$MAIN Syntax: C++ Multiple allowed: no This code block has two main purposes: • Derive new algebraic relationships between parameters, random, effects and other derived variables • Set the initial conditions for model compartments For users who are familiar with NONMEM, $MAIN is similar to $PK. $MAIN is wrapped into a C++ function and compiled / loaded by mrgsolve.
The MAIN function gets called just prior to advancing the system from the current time to the next time for each record in the data set. $MAIN also gets called several times before starting the problem (NEWIND == 0) and just prior to simulating each individual (NEWIND == 1). Finally, $MAIN gets called every time the model initial conditions are queried with init().

This is an alias for $MAIN. ### 2.2.9$ODE

Syntax: C++
Multiple allowed: no
Options: @param

Use $ODE to define model differential equations. For all compartments assign the value of the differential equation to dxdt_CMT where CMT is the name of the compartment. The dxdt_ equation may be a function of model parameters (via $PARAM), the current value of any compartment (CMT) or any user-derived variable.

For example:

[ CMT ] GUT CENT

[ ODE ]
dxdt_GUT = -KA*GUT;
dxdt_CENT = KA*GUT - KE*CENT;

It is important to make sure that there is a dxdt_ expression defined for every compartment listed in $CMT or $INIT, even if it is dxdt_CMT = 0;

The $ODE function is called repeatedly during a simulation run. So it is wise to do as many calculations as possible outside of $ODE, usually in $MAIN. But remember that any calculation that depends on an amount in a compartment and helps determine the dxdt_ expression in a model must be written in $ODE.

New variables may be declared in $ODE. See Section 2.5 for details. For example: $CMT CENT RESP
$PARAM VC = 100, KE = 0.2, KOUT = 2, KIN = 100$ODE
double CP = CENT/VC;
double INH = CP/(IMAX+CP)

dxdt_CENT = -KE*CENT;
dxdt_RESP =  KIN*(1 - INH) - RESP*KOUT;

If the model needs to refer to the current time, use the SOLVERTIME variable.

Notes:

• $ODE is written in C++ syntax; every line must end in ; • There may be only one $ODE block in a model

### 2.2.11$TABLE Syntax: C++ Multiple allowed: no Use $TABLE to interact with parameters, compartment values, and other user-defined variables after the system advances to the next time.

For example:

[ TABLE ]
double CP = CENT/VC;

NOTE mrgsolve formerly had a table() macro for inserting derived values into simulated output. This macro has been deprecated. The only way to insert derived values into the simulated output is via $CAPTURE. NOTE When variables are marked for capture (see Section 2.2.15), the values of those variables are saved at the end of the $TABLE function. This process is carried out automatically by mrgsolve and therefore requires no user intervention.

### 2.2.13$PREAMBLE Syntax: C++ Multiple allowed: no This is the fourth C++ code block. It is called once in two different settings (Chapter 7): 1. Immediately prior to starting the simulation run 2. Immediately prior to calling $MAIN when calculating initial conditions

$PREAMBLE is a function that allows you to set up your C++ environment. It is only called one time during the simulation run (right at the start). The code in this block is typically used to configure or initialize C++ variables or data structures that were declared in $GLOBAL.

For example:

[ PLUGIN ] Rcpp

[ GLOBAL ]
namespace{
Rcpp::NumericVector x;
}

[ PREAMBLE ]
x.push_back(1);
x.push_back(2);
x.push_back(3);

[ MAIN ]
<some code that uses x vector>

In this example, we want to use a numeric vector x and declare it in $GLOBAL so that we can use it anywhere else in the code (the declaration is also made in an unnamed namespace to ensure that the variable is local to the model file). Then, in $PREAMBLE, we put 3 numbers into the vector and we use x in $MAIN. Since $MAIN, $TABLE and (especially) $ODE are called repeatedly as the simulation run progresses, we put the initialization of x in $PREAMBLE to make sure the initialization of x only happens once. Notes: • $PREAMBLE is written in C++ syntax; every line must end in ;
• There may be only one $PREAMBLE block in a model • Like $MAIN, $ODE and $TABLE, double, int and bool variables initialized in $PREAMBLE are actually initialized for global (within the model file) See also: Chapter 7 and Section 2.2.21. ### 2.2.14$PRED

Syntax: C++
Multiple allowed: no

Use $PRED to write a model without differential equations. In this block, write all algebraic expressions for derived parameters, the response, and any other derived output quantities. For example: [ PARAM ] TVE0 = 100, AUC50 = 100, IMAX = 40, AUC = 0 [ PRED ] double E0 = EVE0*exp(ETA(1)); double RESP = E0 - IMAX*AUC/(AUC50+AUC); In this example, the entire model is written in the $PRED block. It is an error to include the following blocks when $PRED is being used: $MAIN, $TABLE, $PKMODEL, $ODE, $CMT, $INIT. See Section 3.7 for additional information regarding data sets in use with $PRED block.

### 2.2.15$CAPTURE Syntax: text Multiple allowed: yes Options: @annotated, @etas This is a block to identify variables that should be captured in the simulated output. The @etas option is new with version 1.0.8. For example: [ PARAM ] A = 1, B = 2 [ MAIN ] double C = 3; bool yes = true; [ CAPTURE ] A B C yes This construct will result in four additional columns in the simulated output with names A, B, C, and yes. Users can also rename captured variables by providing a newname = oldname specification. $PARAM WT = 70, THETA1 = 2.2

$MAIN double CL = THETA1*pow(WT/70,0.75)*exp(ETA(1));$OMEGA 1

$CAPTURE WEIGHT = WT TVCL = THETA2 CL ETA(1) In this example, the names of the captured data items will be WEIGHT,TVCL,CL,ETA_1. Users can use the capture type to declare variables in $MAIN and $TABLE. capture types are really doubles, but using that type will signal mrgsolve to automatically capture that value. For example: $PARAM VC = 300

$CMT CENT$TABLE
capture DV = (CENT/VC);

Since we used type capture for DV, DV will show up as a column in the simulated data.

Annotated example:

$MAIN double CLi = TVCL*exp(ECL);$TABLE
double DV = (CENT/VC)*exp(PROP);

$CAPTURE @annotated CLi : Individual clearance (L/hr) DV : Plasma concentration (mcg/ml) Use the @etas option to specify an expression that evaluates to integers which identify the ETA number(s) to capture. For example $OMEGA 0.1 0.2 0.3

$CAPTURE @etas 1:LAST CL DV Will capture all three ETAs into the simulation outputs. NONMEM-style names will be used (e.g. ETA1). This is in contrast to the names that are generated when ETA(1) is captured; in this case, the parens are removed and replaced with _ so the output name is ETA_1. mrgsolve will provide LAST (and last) which will evaluate to the maximum number of ETAs in the problem (in this case, 3). You can write any valid expression in @etas. When the expression is evaluated, it will be evaluated in the model environment created through the $ENV block (Section 2.2.26).

Tip

In addition to listing model variables for otput in the $CAPTURE block, users can “dynamically” capture model variables through the capture argument to mread(). These may be model parameters (listed in $PARAM) or user-defined model variables.

To get a listing of available user-defined variables, you can use

mod <- house()
as.list(mod)$cpp_variables ### 2.2.16$OMEGA

Syntax: text
Multiple allowed: yes
Options: @annotated, @block, @correlation, @labels, @name, @object, @as_object

See ?modMATRIX for more details about options for this block.

Use this block to enter variance/covariance matrices for subject-level random effects drawn from multivariate normal distribution. All random effects are assumed to have mean of 0. Off diagonal elements for block matrices are assumed to be correlation coefficients if the @correlation option is used (see below).

By default, a diagonal matrix is assumed. So:

$OMEGA 1 2 3 will generate a 3x3 omega matrix. A block matrix may be entered by using block = TRUE. So: $OMEGA @block
0.1 0.02 0.3

will generate a 2x2 matrix with covariance 0.02.

A 2x2 matrix where the off-diagonal element is a correlation, not a covariance can be specified like this:

$OMEGA @correlation 0.1 0.67 0.3 Here, the correlation is 0.67. mrgsolve will calculate the covariances and substitute these values. The matrix will be stored and used with these covariances, not the correlation. A name can be assigned to each matrix: $OMEGA @name PK @block
0.2 0.02 0.3

$OMEGA @name PD 0.1 0.2 0.3 0.5 to distinguish between multiple $OMEGA blocks and to facilitate updating later. The model in the preceding example will have two $OMEGA matrices: 2x2 and 4x4. Labels can be assigned which function as aliases for the different ETAs $OMEGA @block @labels ETA_CL ETA_V
0.1 0.05 0.2

The number of labels should match the number of rows (or columns) in the matrix.

Annotated example (diagonal matrix):

$OMEGA @annotated ECL: 0.09 : ETA on clearance EVC: 0.19 : ETA on volume EKA: 0.45 : ETA on absorption rate constant Annotated example (block matrix): $OMEGA @annotated @block
ECL: 0.09 : ETA on clearance
EVC: 0.001 0.19 : ETA on volume
EKA: 0.001 0.001 0.45 : ETA on absorption rate constant

The @object and @as_object options should name or return a square numeric matrix. If rownames are included in the matrix, then they will be used to form labels for the realized ETAs. For example, we can initialize a very large $OMEGA matrix with $OMEGA @as_object
n <- 20
lbl <- paste0("ETA_", LETTERS[1:n])
matrix(0, 20, 20, dimnames = list(lbl, lbl))

Note: this only initializes the matrix; you will (likely) need to update it with meaningful values after the model is loaded (Chapter 5).

### 2.2.18$SET Syntax: R Multiple allowed: no Use this code block to set different options for the simulation. Use a name = value format, where value is evaluated by the R interpreter. Most of the options that can be entered in $SET are passed to update.

For example:

[ SET ] end = 240, delta = 0.5, req = "RESP"

Here, we set the simulation end time to 240, set the time difference between two adjacent time points to 0.25 time units, and request only the RESP compartment in the simulated output.

### 2.2.19$GLOBAL Syntax: C++ Multiple allowed: no The $GLOBAL block is for writing C++ code that is outside of $MAIN, $ODE, and $TABLE. There are no artificial limit on what sort of C++ code can go in $GLOBAL.
However there are two more-common uses:

1. Write #define preprocessor statements
2. Define global variables, usually variables other than double, bool, int (see Section 2.5)

Preprocessor directives Preprocessor #define directives are direct substitutions that the C++ preprocessor makes prior to compiling your code.

For example:

[ GLOBAL ]
#define CP (CENT/VC)

When this preprocessor directive is included, everywhere the preprocessor finds a CP token it will substitute (CENT/VC). Both CENT and VC must be defined and the ratio of CENT to VC will be calculated depending on whatever the current values are. Notice that we included parentheses around (CENT/VC).
This makes sure the ratio between the two is taken first, before any other operations involving CP.

Declaring global variables Sometimes, you may wish to use global variables and have more control over how they get declared.

$GLOBAL bool cure = false; With this construct, the boolean variable cure is declared and defined right as the model is compiled. Declare in bulk If you have a large number of variables to declare, you can do that in bulk in $GLOBAL. For example, we know we will need a long list of double precision variable for covariate modeling, we can declare them all at once like this:

[ global ]
double TVCL, TVV2, TVQ = 0, TVV3 = 0;

[ main ]

TVCL = THETA1 * pow(WT / 70.0, 0.75);
TVV2 = THETA2 * WT / 70.0;
TVQ  = THETA3 * pow(WT / 70.0, 0.75);
TVV3 = THETA4 * WT / 70.0;

This isn’t a terribly long list, but let’s pretend it is to illustrate how to do this. You can also declare int, bool and other types like this. I’ve initialized the last two (TVVQ and TVV3) to illustrate how to do this. It’s good practice to do that but not necessary as long as everything gets initialized to something before they are used.

$PKMODEL cmt="GUT CENT PERIPH", depot = TRUE Specifying three compartments with depot=TRUE implies ncmt=2. Notice that a separate $CMT block is not appropriate when cmt is specified in $PKMODEL. ### 2.2.21$PLUGIN

Syntax: text
Multiple allowed: no

Plugins are a way to add extensions to your mrgsolve model. Plugins can either link your model to external libraries (like boost or Rcpp) or they can open up access to additional functionality provided by mrgsolve itself.

Plugins are listed and discussed in more detail in Chapter 9.

Usage

To invoke a plugin, list the plugin name in the code block. For requesting Rcpp headers in your model, call

$PLUGIN Rcpp Available plugins The following plugins make additional mrgsolve-specific functionality available • mrgx: extra C++ functions (see below) • tad: track time after dose in your model • N_CMT: get the number of a compartment The following plugins give you a different model specification experience • autodec: mrgsolve will find assignments in your model and automatically declare them as doubles • nm-vars: a NONMEM look and feel for compartmental models; use F1, A(1) and DADT(1) rather than F_GUT, GUT and dxdt_GUT The following plugin lets you customize how your model is compiled • CXX11: compile your model with C++11 standard; this adds the compiler flag -std=c++11 The following plugins will let you link to external libraries • Rcpp: include Rcpp headers int your model • BH: include boost headers in your model • RcppArmadillo: include Armadillo headers in your model Note that Rcpp, RcppArmadillo and BH only allow you to link to those headers. To take advantage of that, you will need to know how to use Rcpp, boost etc. For the BH plugin, no headers are included for you; you must include the proper headers you want to use in $GLOBAL.

mrgx This is a general collection of functions that we made available.

Functions provided by mrgx:

• T get<T>(std::string <pkgname>, std::string <objectname>)
• This gets an object of any Rcpp-representable type (T) from any package
• T get<T>(std::string <objectname)
• This gets an object of any Rcpp-representable type (T) from .GlobalEnv
• T get<T>(std::string <objectname>, databox& self)
• This gets an object of any Rcpp-representable type (T) from $ENV • double rnorm(double mean, double sd, double min, double max) • Simulate one variate from a normal distribution that is between min and max • double rlognorm(double mean, double sd, double min, double max) • Same as mrgx::rnorm, but the simulated value is passed to exp after simulating • Rcpp::Function mt_fun() • Returns mrgsolve::mt_fun; this is usually used when declaring a R function in $GLOBAL
• Example: Rcpp::Function print = mrgx::mt_fun();

IMPORTANT All of these functions are in the mrgx namespace. So, in order to call these functions you must include mrgx:: namespace identifier to the front of the function name. For example, don’t use rnorm(50,20,40,140); use mrgx::rnorm(50,20,40,140).

#### 2.2.21.1 Some examples

Get a numeric vector from $ENV [ PLUGIN ] Rcpp mrgx [ ENV ] x <- c(1,2,3,4,5) [ GLOBAL ] Rcpp::NumericVector x; [ PREAMBLE ] x = mrgx::get<Rcpp::NumericVector>("x", self); Get the print function from package:base $PLUGIN Rcpp mrgx

$GLOBAL Rcpp::Function print = mrgx::mt_fun();$PREAMBLE
print = mrgx::get<Rcpp::Function>("base", "print");

$MAIN print(self.rown); Note that we declare the print in $GLOBAL and use the mt_fun() place holder.

Simulate truncated normal variables This simulates a weight that has mean 80, standard deviation 20 and is greater than 40 and less than 140.

$PLUGIN Rcpp mrgx$MAIN
if(NEWIND <=1) {
double WT = mrgx::rnorm(80,20,40,140);
}

### 2.2.22$THETA Syntax: text Multiple allowed: yes Options: @annotated, @name, @object, @as_object Use this code block as an efficient way to add to the parameter list where names are determined by a prefix and a number. By default, the prefix is THETA and the number sequentially numbers the input values. For example: [ THETA ] 0.1 0.2 0.3 is equivalent to $PARAM THETA1 = 0.1, THETA2 = 0.2, THETA3 = 0.3

Annotated example:

$THETA @annotated 0.1 : Typical value of clearance (L/hr) 0.2 : Typical value of volume (L) 0.3 : Typical value of ka (1/hr) To change the prefix, use @name option $THETA @name theta
0.1 0.2 0.3

would be equivalent to

[ PARAM ] theta1 = 0.1, theta2 = 0.2, theta3 = 0.3

The @object and @as_object options should name or return an unnamed vector of parameter values (names are ignored).

### 2.2.23$NMXML Syntax: R Multiple allowed: yes The $NMXML block lets you read and incorporate results from a NONMEM run into your mrgsolve model. For example

$NMXML run = 1001 project = "../model/nonmem" root = "cppfile" From the NONMEM run, THETA will be imported into your parameter list (see Section 2.2.4 and Section 1.1), OMEGA will be captured as an $OMEGA block (Section 2.2.16) and SIGMA will be captured as a $SIGMA block (Section 2.2.17). Users may optionally omit any one of these from being imported. $NMXML contains a project argument and a run argument. By default, the estimates are read from from the file project/run/run.xml. That is, it is assumed that there is a directory named run that is inside the project directory where $NMXML will find run.xml. Your NONMEM run directories may not be organized in a way that is compatible with this default. In that case, you will need to provide the file argument, which should be the path to the run.xml file, either as a full path or as a path relative to the current working directory. Once the model object is obtained, the path to the xml file that formed the source for imported parameters can be retrieved by coercing the model object to list and looking for nm_import: mod <- modlib("1005", compile = FALSE) as.list(mod)$nm_import

For help on the arguments / options for $NMXML, please see the ?nmxml help topic in your R session after loading the mrgsolve package. An example There is a NONMEM run embedded in the mrgsolve package path <- file.path(path.package("mrgsolve"),"nonmem") list.files(path, recursive=TRUE) . [1] "1005/1005.cat" "1005/1005.coi" "1005/1005.cor" "1005/1005.cov" . [5] "1005/1005.cpu" "1005/1005.ctl" "1005/1005.ext" "1005/1005.grd" . [9] "1005/1005.lst" "1005/1005.phi" "1005/1005.shk" "1005/1005.shm" . [13] "1005/1005.tab" "1005/1005.xml" "1005/1005par.tab" "1005/INTER" . [17] "2005/2005.ext" We can create a mrgsolve control stream that will import THETA, OMEGA and SIGMA from that run using the $NMXML code block.

// inline/nmxml-1.cpp

$NMXML run = 1005 project = path root = "cppfile" olabels = c("ECL", "EVC", "EKA") slabels = c("PROP", "ADD")$MAIN
double CL = THETA1*exp(ECL);
double V2 = THETA2*exp(EVC);
double KA = THETA3*exp(EKA);
double Q = THETA4;
double V3 = THETA5;

$PKMODEL ncmt=2, depot=TRUE$CMT GUT CENT PERIPH

$TABLE double CP = (CENT/V2)*(1+PROP) + ADD/5;$CAPTURE CP

$SET delta=4, end=96 NOTE: in order to use this code, we need to install the xml2 package. mod <- mread("inline/nmxml-1.cpp", compile = FALSE) mod . . . --------------- source: nmxml-1.cpp --------------- . . project: /Users/kyleb/git...-guide/inline . shared object: nmxml-1.cpp-so-1058b51a4e6d6 <not loaded> . . time: start: 0 end: 96 delta: 4 . add: <none> . . compartments: GUT CENT PERIPH [3] . parameters: THETA1 THETA2 THETA3 THETA4 THETA5 . THETA6 THETA7 [7] . captures: CP [1] . omega: 3x3 . sigma: 2x2 . . solver: atol: 1e-08 rtol: 1e-08 maxsteps: 20k . ------------------------------------------------------ param(mod) . . Model parameters (N=7): . name value . name value . THETA1 9.51 | THETA5 113 . THETA2 22.8 | THETA6 1.02 . THETA3 0.0714 | THETA7 1.19 . THETA4 3.47 | . . revar(mod) .$omega
. $... . [,1] [,2] [,3] . ECL: 0.21387884 0.12077020 -0.01162777 . EVC: 0.12077020 0.09451047 -0.03720637 . EKA: -0.01162777 -0.03720637 0.04656315 . . .$sigma
. $... . [,1] [,2] . PROP: 0.04917071 0.0000000 . ADD: 0.00000000 0.2017688 root argument: please use the root = "cppfile" argument going forward. As of mrgsolve 0.11.0, we added an argument called root to $NMXML that tells mrgsolve the location where it should read the xml file from. The default behavior is the "working" directory. When this is the case, mrgsolve assumes that the the xml file can be found relative to the “working” directory. The only other value that root can take is "cppfile". When root = "cppfile", then mrgsolve will look for the xml file in a directory that is relative to where the model source code file is located. Please take a look at Section 2.2.24 for more discussion and examples.

Starting with version 1.0.8, when the NONMEM run is found using the run + project arguments, you can pass run = "@cppstem" and mrgsolve will look for the NONMEM run number that matches the stem of the current mrgsolve file. For example, if you are coding a model to match NONMEM run 101.ctl, then call your mrgsolve model file 101.cpp or 101.mod; just match the file stem. Then write

101.cpp
$NMXML run = "@cppstem" project = <nonmem-run-directory> mrgsolve will check the stem of the mrgsolve model file and look for the NONMEM .xml file in <nonmem-run-directory>/101/101.xml Please see the ?nmxml help topic for more information on arguments that can be passed to $NMXML.

### 2.2.24$NMEXT Syntax: R Multiple allowed: yes Like $NMXML, $NMEXT allows the import of $THETA, $OMEGA, and $SIGMA from your NONMEM run into your mrgsolve model, but the estimates are read from the .ext file output. $NMEXT is able to import the NONMEM estimates much faster than $NMXML when loading from sampling based methods (mainly METHOD=BAYES).

We can load the same mode as the $NMXML example above with the following code: // inline/nmext-1.cpp$NMEXT
run = 1005
project = file.path(path.package("mrgsolve"),"nonmem")
root = "cppfile"

olabels = c("ECL", "EVC", "EKA")

$MAIN double CL = THETA1*exp(ECL); double V2 = THETA2*exp(EVC); double KA = THETA3*exp(EKA); double Q = THETA4; double V3 = THETA5;$PKMODEL ncmt=2, depot=TRUE

$CMT GUT CENT PERIPH$TABLE
double CP = (CENT/V2)*(1+PROP) + ADD/5;

$CAPTURE CP$SET delta = 4, end = 96
mod <- mread("inline/nmext-1.cpp", quiet=TRUE)

param(mod)
.
.  Model parameters (N=7):
.  name   value  . name   value
.  THETA1 9.51   | THETA5 113
.  THETA2 22.8   | THETA6 1.02
.  THETA3 0.0714 | THETA7 1.19
.  THETA4 3.47   | .      .

Once the model object is obtained, the path to the ext file that formed the source for imported parameters can be retrieved by coercing the model object to list and looking for nm_import. See the $NMXML topic in Section 2.2.23 for an example of equivalent functionality. IMPORTANT: while $NMEXT works very similarly to $NMXML, there is one key difference between the two: when using $NMEXT, you will always get one $OMEGA matrix and one $SIGMA matrix, regardless of the block structure used in the NONMEM control stream. When using $NMXML, you get the same block structure in the mrgsolve model as the NONMEM model. For example: in the NONMEM control stream, $OMEGA is one 2x2 matrix and one 3x3 matrix. Importing those estimates with $NMEXT will give you one 5x5 matrix, while importing the estimates with $NMXML will give you a list of one 2x2 and one 3x3 matrix.

root argument: please use the root = "cppfile" argument going forward.

As of mrgsolve 0.11.0, we added an argument called root to $NMEXT that tells mrgsolve the location where it should read the ext file from. The default behavior is the "working" directory. When this is the case, mrgsolve assumes that the the ext file can be found relative to the “working” directory. The only other value that root can take is "cppfile". When root = "cppfile", then mrgsolve will look for the ext file in a directory that is relative to where the model source code file is located. We recommend that users start using the root argument and set it to “cppfile”. This will eventually become the default. For example: $NMEXT
run = 1005
project = "../../model/pk"
root = "cppfile"

This tells mrgsolve to find the nonmem run back two directories (../../) and then into model --> pk --> 1005 relative to where the mrgsolve model file is located. This is in contrast to the previous expected behavior that the path should be relative to the current working directory.

We also note here that when users pass an absolute path, the relative path doesn’t matter at all. Users are free to use pathing tools to generate the absolute path to the project. For example, the here::here() function can be used like this

$NMEXT run = 1005 project = here::here("model/pk") When used in the proper context, here() will generate an absolute path from your projects root directory to the nonmem project directory. Please refer to the here::here() documentation for proper use of this function. Starting with version 1.0.8, when the NONMEM run is found using the run + project arguments, you can pass run = "@cppstem" and mrgsolve will look for the NONMEM run number that matches the stem of the current mrgsolve file. See details and example in Section 2.2.23. See the ?nmext R help topic for arguments that can be passed to $NMEXT block. Notably, the user can select the function to read in the ext file.By default, mrgsolve will try to load data.table and use the fread function. If data.table can’t be loaded, then mrgsolve will use utils::read.table.

### 2.2.25$INCLUDE Syntax: text Multiple allowed: no To include your own header file(s) in a model use $INCLUDE

$INCLUDE mystuff.h otherstuff.h or $INCLUDE
mystuff.h,  otherstuff.h

mrgsolve will insert proper #include preprocessor directives into the C++ code that gets compiled.

Requirements

• All header files listed in $INCLUDE are assumed (expected) to be under the project directory; don’t use $INCLUDE for header files that are in any other location
• An error is generated if the header file does not exist
• An error is generated if any quotation marks are found in the file name (don’t use quotes around the file name; mrgsolve will take care of that)
• A warning is issued if the header file does not end in .h
• When the header file is changed (MD5 checksum changes), the model will be forced to be rebuilt (recompiled) when mread() or mcode() (but not mread_cache() or mcode_cache()) is called; this feature is only available for header files listed in $INCLUDE (see below) • Do not use $INCLUDE to include Rcpp, boost, RcppArmadillo or RcppEigen headers; use the appropriate $PLUGIN instead For applications that don’t fit into the requirements listed above, users can always include header files in the model in $GLOBAL like this:

$GLOBAL #include "/Users/me/libs/mystuff.h" But be careful when doing this: if there are changes to mystuff.h but not to any other part of the model specification, the model may not be fully compiled when calling mread. In this case, always use preclean=TRUE argument to mread to force the model to be built when calling mread. ### 2.2.26$ENV

Syntax: R
Multiple allowed: no

This block is all R code (just as you would code in a stand-alone R script. The code is parsed and evaluated into a new environment when the model is compiled.

For example:

$ENV Sigma <- cmat(1,0.6,2) mu <- c(2,4) cama <- function(mod) { mod %>% ev(amt=100, ii=12, addl=10) %>% mrgsim(obsonly=TRUE,end=120) } Objects inside $ENV can be utilized in different C++ functions (see Section 2.2.21) or other parts of the simulation process.

Objects inside $ENV will also be made available when parsing different parts of the model specification file. For example, we can initialize an $OMEGA block with the Sigma object in the $ENV block example above: $OMEGA @object Sigma

Another example is providing constants or functions to simplify code in the $PARAM block $ENV
rescale <- 1/1000

convert <- function(x) x * 5 / 166.2

$PARAM p = 500 * rescale q = convert(2.51) You can also make your model specification file depend on an R package; for example, if you want to use the here package to render a path to your NONMEM model $ENV
if(!requireNamespace("here")) {
stop("the here package is required to load this model.")
}

$MAIN if(CMT==2) { // .... } In the example, GUT, CENT and PERIPH are the amounts in the respective compartments and CMT refers to the value of CMT in the data record / data set. CMT is an alias for self.cmt. ### 2.3.7AMT AMT is the current value of dose amount. AMT is an alias for self.amt. ### 2.3.8NEWIND NEWIND is a new individual indicator, taking the following values: • 0 for the first event record of the data set • 1 for the first event record of a subsequent individual • 2 for subsequent event record for an individual For example: [ GLOBAL ] int counter = 0; [ MAIN ] if(NEWIND <=1) { counter = 0; } NEWIND is an alias for self.newind. ### 2.3.9SS_ADVANCE This is a bool data item (either true or false) which is always false unless mrgsolve is currently advancing the system to steady state (Chapter 8); then it will be true. This variable is only available in the $ODE (or $DES) block. Use this variable to modify the calculation of your model differential equations while the system is advancing to steady state. One use case is when you have an accumulation compartment to calculate AUC but you want to halt accumulation while the system is working toward steady state $ODE
dxdt_CENT = -(CL/V) * CENT;
dxdt_AUC = CENT/V;
if(SS_ADVANCE) dxdt_AUC = 0;

### 2.3.10simeta()

The simeta() function can be used to re-simulate ETA values. For example,

$MAIN simeta(); will re-simulate all ETA(n). ### 2.3.11simeps() simeps() works like simeta(), but all EPS(n) values are re-simulated rather than the ETA values. ### 2.3.12self self is an object that gets passed to your model function that contains data items and functions that can be called. It is a struct (see the source code here). A partial list of members are documented in the following sections. self functions include • self.tad() • self.mtime() • self.mevent() • self.stop() • self.stop_id() • self.stop_id_cf() self data members include • self.id (ID) • self.amt (AMT) • self.cmt (CMT) • self.evid (EVID) • self.time (TIME) • self.newind (NEWIND) • self.nid • self.idn • self.rown • self.nrow • self.envir (see note below) This information is provided for transparency and is not exhaustive. We provide an interface to these data items through pre-processor directives with simpler names (e.g. EVID will translate to self.evid) and users are encouraged to use the simpler, API name. ### 2.3.13self.time The current data set TIME. ### 2.3.14self.cmt The current compartment number regardless of whether it was given as cmt or CMT in the data set. There is no alias for self.cmt. For example: $TABLE

double DV = CENT/VC + EPS(1);

if(self.cmt==3) DV = RESPOSE + EPS(2);

### 2.3.15self.amt

The current amt value regardless of whether it was given as amt or AMT in the data set. There is no alias for self.amt.

[ PREAMBLE ]
double last_dose = 0;

[ MAIN ]

if(EVID==1) {
last_dose = self.amt;
}

### 2.3.16self.nid

The number of IDs in the data set.

### 2.3.17self.idn

The current ID number. Numbers start at 0 and increase by one to self.nid-1. So if you want to test for the last ID in the output data set, you would write:

[ table ]

if(self.idn == (self.nid-1)) {
// do something ....
}

### 2.3.18self.nrow

The total number of rows in the output data set.

### 2.3.19self.rown

The current row number. Numbers start at 0 and increase by one to self.rown-1. So if you want to test for the last row in the output data set, you would write:

[ table ]

if(self.rown == (self.nrow-1)) {
// do something ....
}

### 2.3.20self.envir

This item is a null pointer that can be cast to an Rcpp::Environment only when the Rcpp plugin is invoked (see Section 9.5).

For example

[ plugin ] Rcpp mrgx

[ preamble ]
Rcpp::Environment env = mrgx::get_envir(self);

### 2.3.21self.tad()

This is a function that calculates and returns the time after the most recent dose event record (any record with EVID equal to 1 or 4). self.tad() will return -1 when it is called before the first dose record within an individual (NEWIND <= 1).

This function should be called in $MAIN every time $MAIN is called.

For example:

$MAIN double TAD = self.tad(); Do not make this calculation depend on any test or condition; it must be called every time $MAIN is called in order to see every dose.

### 2.3.22self.mtime(<time>)

This is a function that creates a modeled even time. The only argument is a numeric time in the future for the current individual that indicates when a discontinuity should be included in the simulation.

When when self.mtime() is called, a new record is added to the current individual’s record set with EVID set to 2. This means that when the system advances to this record (this time), the differential equation solver will reset and restart, creating the discontinuity. The function returns the time of this event so the user can work with it in subsequent code.

For example,

$PARAM change_point = 5.13;$MAIN

double KA = 1.1;

double mt1 = self.mtime(change_point);

if(TIME >= mt1) KA = 2.1;  

### 2.3.23self.mevent(<time>, <evid>)

Related to self.mtime() (Section 2.3.22), except you can set a specific EVID for this intervention and you track the change time via EVID rather than the time. You’d use this if you want to anchor several events in the future and be sure you can distinguish between them. For example

[ main ]
self.mevent(change_point1, 33);
self.mevent(change_point2, 34);

Now, you can look for EVID==33 to know you have hit change_point1 or EVID==34 to indicate that you have hit change_point2. Notice that self.mevent() doesn’t return any value. You will know that you’ve hit the change time by checking EVID.

This self function is available to be called from $PREAMBLE, $MAIN, and $TABLE. When this function is called, the entire problem is stopped upon processing the next simulation record. This might be called when something really bad happened and you just want to stop the simulation with an error. ### 2.3.25self.stop_id() This self function is available to be called from $PREAMBLE, $MAIN, and $TABLE. When this function is called, processing of the current individual is stopped and missing values (NA_real) are filled in for remaining compartment’ and capture outputs.

This might be called when some condition is reached in the current individual that indicates either that the rest of the outputs are inconsequential or there was a problem with this particular individual.

See also self.stop_id() and self.stop().

1 2 3

$TABLE double STD = sqrt(SIGMA(1) + pow(F,2)*SIGMA(2)); ### 2.3.31table(<name>) This macro has been deprecated. Users should not use code like this: [ TABLE ] table(CP) = CENT/VC; But rather this: $TABLE
double CP = CENT/VC;

$CAPTURE CP See: Section 2.2.11 and also Section 2.2.15. ### 2.3.32F_CMT For the CMT compartment, sets the bioavailability fraction for that compartment. Example: $MAIN
F_CENT = 0.7;

### 2.3.33ALAG_CMT

For the CMT compartment, sets the lag time for doses into that compartment.

Example:

$MAIN ALAG_GUT = 0.25; ### 2.3.34R_CMT For the CMT compartment, sets the infusion rate for that compartment. The infusion rate is only set via R_CMT when rate in the data set or event object is set to -1. Example: $MAIN
R_CENT = 100;

### 2.3.35Rn

Only available when the nm-vars plugin is invoked; see R_CMT.

### 2.3.36D_CMT

For the CMT compartment, sets the infusion duration for that compartment.
The infusion duration is only set via D_CMT when rate in the data set or event object is set to -2.

Example:

$MAIN D_CENT = 2; ### 2.3.37 NONMEM-like syntax There is a NONMEM-like syntax that lets you write F1 rather than F_GUT, D2 rather than D_CENT as well as other variables that are commonly used in a NONMEM control stream. To make this syntax available, the nm-vars plugin must be invoked. See Section 9.2 on the nm-vars plugin for details. ## 2.4 Reserved words Reserved words cannot be used as names for parameters, compartments or other derived variables in the model. Note that some of these words are “reserved” for you to use in your data set. ID amt cmt ii ss evid addl rate time SOLVERTIME table ETA EPS AMT CMT ID TIME EVID simeps self simeta NEWIND DONE CFONSTOP DXDTZERO CFONSTOP INITSOLV _F _R _ALAG SETINIT report _VARS_ VARS SS_ADVANCE AMT CMT II SS ADDL RATE THETA pred_CL pred_VC pred_V pred_V2 pred_KA pred_Q pred_VP pred_V3 double int bool capture Other reserved words depend on the compartment names in your model. For example, if you have a compartment called CENT in the model, then the the following will be reserved F_CENT R_CENT D_CENT ALAG_CENT N_CENT ## 2.5 Derive new variables New C++ variables may be derived in $GLOBAL, $PREAMBLE $MAIN, $ODE and $TABLE. Because these are C++ variables, the type of variable being used must be declared. For the vast majority of applications, the double type is used (double-precision numeric value).

$MAIN double CLi = TVCL*exp(ETA(1)); We want CLi to be a numeric value, so we use double. To initialize a boolean variable (true / false), write $MAIN
bool cure = false;

### 2.5.1 Special handling for double, int, bool

When variables of the type double, int, and bool are declared and initialized in $PREAMBLE, $MAIN, $ODE, $TABLE, mrgsolve will detect those declarations, and modify the code so that the variables are actually declared once in $GLOBAL not in $MAIN, $ODE, or $TABLE. This is done so that variables declared in one code block (e.g. $MAIN) can be read and modified in another code block (e.g. $TABLE).

For example, in the following code:

$MAIN double CLi = TVCL*exp(ETA(1)); a double-precision numeric variable is created (CLi) in the $MAIN block. When mrgsolve parses the model file, this code gets translated to

$GLOBAL namespace { double CLi; }$MAIN
CLi = TVCL*exp(ETA(1));

That is, CLi is declared in $GLOBAL in an unnamed namespace so that variables like this are global variables within the model file only. This way, we can still read the CLi variable in $TABLE:

$MAIN double CLi = TVCL*exp(ETA(1)); double VCi = TVVC*exp(ETA(2));$TABLE
double KEi = CLi/VCi;

$CAPTURE KEi To declare a variable that is local to a particular code block: $MAIN

localdouble CLi = TVCL*exp(ETA(1));

The localdouble type is still just a double-precision variable. The difference is that it is protected from this re-declaration process and the variable will be local to (in this case) the $MAIN block. ### 2.5.2 Using other types globally As we noted in the previous section, double, int, and bool are processed in a special way so that they are by default global to the file. Many times we want to work with other variable types in a global manner. Whenever you want a data structure to be accessible across functions (e.g. $MAIN, $TABLE, etc.) they should be declared in $GLOBAL, optionally in an unnamed namespace.

For example:

[ GLOBAL ]
std::vector<double> myvec;

or

[ GLOBAL ]
namespace {
std::vector<double> myvec;
}

In case that object needs some configuration prior to starting the problem, use $PREAMBLE to do that work [ GLOBAL ] std::vector<double> myvec; [ PREAMBLE ] myvec.assign(3,1); ## 2.6 Random number generation Users can simulate random numbers inside a model file using functions that are similar to the functions you’d normally use in R (e.g. rnorm() and runif()). This functionality is provided by Rcpp and therefore requires using the Rcpp plugin (see Chapter 9 and Section 9.5). Rcpp provides these functions inside the R namespace so you will have to prefix the function call with R::. As an example, to make a draw from Uniform (0,1) [ plugin ] Rcpp [ error ] double draw = R::runif(0,1); Note that the 0 gets used as min and the 1 gets used as max; we didn’t pass n here and draw is a single number (not a vector like you’d get from (runif(100, 0, 1) on your R console). So in general, these functions work like their R counterparts, but without the n argument. Another example showing how to draw from binomial distribution with probability 0.5 [ plugin ] Rcpp [ error ] double draw = R::rbinom(1, 0.5); Here, the 1 is used as size (not n) and 0.5 is used as prob. Other helpful functions could be R::rnorm() or R::rlnorm() but you can call any of the r functions as well as the corresponding dpq functions through this R namespace. Documentation, including functions to call and arguments can be found in the Rcpp API docs http://dirk.eddelbuettel.com/code/rcpp/html/namespaceR.html ## 2.7 Examples The following sections show example model specification. The intention is to show how the different blocks, macros and variables can work together to make a functional model. Some models are given purely for illustrative purpose and may not be particularly useful in application. ### 2.7.1 Simple PK model Notes: • Basic PK parameters are declared in $PARAM; every parameter needs to be assigned a value
• Two compartments GUT and CENT are declared in $CMT; using $CMT assumes that both compartments start with 0 mass
• Because we declared GUT and CENT as compartments, we write dxdt_ equations for both in $ODE • In $ODE, we refer to parameters (CL/VC/KA) and the amounts in each compartment at any particular time (GUT and CENT)
• $ODE should be C++ code; each line ends in ; • We derive a variable called CP in $TABLE that has type capture; mrgsolve will enter the CP name into the $CAPTURE block list $PARAM CL = 1, VC = 30, KA = 1.3

$CMT GUT CENT$ODE

dxdt_GUT = -KA*GUT;
dxdt_CENT = KA*GUT - (CL/VC)*CENT;

$TABLE capture CP = CENT/VC; This model can also be written without differential equations [ PARAM ] CL = 1, VC = 30, KA = 1.3 [ PKMODEL ] cmt = "CMT GUT CENT", depot = TRUE$TABLE
capture CP = CENT/VC;

### 2.7.2 PK/PD model

Notes:

• We use a preprocessor #define directive in $GLOBAL; everywhere in the model where a CP token is found, the expression (CENT/VC) … with parentheses … is inserted • We write the initial value for the RESP compartment in $MAIN as a function of two parameters KIN/KOUT
• A new variable - INH- is declared and used in $ODE • Since CP is defined as CENT/VC, we can “capture” that name/value in $CAPTURE
• Both $MAIN and $ODE are C++ code blocks; don’t forget to add the ; at the end of each statement
$PARAM CL = 1, VC = 30, KA = 1.3 KIN = 100, KOUT = 2, IC50 = 2$GLOBAL
#define CP (CENT/VC)

$CMT GUT CENT RESP$MAIN
RESP_0 = KIN/KOUT;

$ODE double INH = CP/(IC50+CP); dxdt_GUT = -KA*GUT; dxdt_CENT = KA*GUT - (CL/VC)*CENT; dxdt_RESP = KIN*(1-INH) - KOUT*RESP;$CAPTURE CP

### 2.7.3 Population PK model with covariates and IOV

Notes:

• Use $SET to set the simulation time grid from 0 to 240 by 0.1 • There are two $OMEGA matrices; we name them IIV and IOV
• The IIV “etas” are labeled as ECL/EVC/EKA; these are aliases to ETA(1)/ETA(2)/ETA(3). The IOV matrix is unlabeled; we must refer to ETA(4)/ETA(5) for this
• Because ETA(1) and ETA(2) are labeled, we can “capture” them as ECL and EVC
• We added zeros for both $OMEGA matrices; all the etas will be zero until we populate those matrices (Section 5.3) $PARAM TVCL = 1.3, TVVC=28, TVKA=0.6, WT=70, OCC=1

$SET delta=0.1, end=240$CMT GUT CENT

$MAIN double IOV = IOV1 if(OCC==2) IOV = IOV2; double CLi = exp(log(TVCL) + 0.75*log(WT/70) + ECL + IOV); double VCi = exp(log(TVVC) + EVC); double KAi = exp(log(TVKA) + EKA);$OMEGA @name IIV @labels ECL EVC EKA
0 0 0
$OMEGA @name IOV @labels IOV1 IOV2 0 0$SIGMA 0

$ODE dxdt_GUT = -KAi*GUT; dxdt_CENT = KAi*GUT - (CLi/VCi)*CENT;$TABLE
capture CP = CENT/VCi;

\$CAPTURE IOV ECL EVC