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:

Code in a separate file and source into your R script
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")

to load that model.

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.

2.1.3 Comments

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
$INPUT
$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 12.1 for an example.

2.2.4 `$PARAM`

Syntax: R
Multiple allowed: yes
Options: @annotated, @covariates, @tag, @input, @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.

The @input option adds a tag to the parameters in the block so that input data sets will be checked for these parameters when check_data_names() is called. Use the @tag option and supply a user-defined tag that can also be used for checking data sets and parameters.

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 MWT = 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.5, Section 2.2.24 and Section 2.2.6.

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

2.2.5 `$INPUT`

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

This block operates just like $PARAM (Section 2.2.4), in that it is a way to introduce parameters in to the model. Additionally, the @input tag is added to the model parameters and input data sets will be checked for these parameters when check_data_names() is called.

2.2.6 `$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.24.

Notes:

Multiple blocks are allowed
Values are evaluated by the R interpreter

2.2.7 `$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.8 `$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().

New variables may be declared in $MAIN. See Section 2.5 for details.

Examples:

[ CMT ] CENT RESP

[ PARAM ] KIN = 100, KOUT = 2, CL = 1, VC = 20

[ MAIN ]

RESP_0 = KIN/KOUT;

double ke = CL/VC;

2.2.9 `$PK`

This is an alias for $MAIN.

2.2.10 `$ODE`

Syntax: C++
Multiple allowed: yes 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 `$DES`

This is an alias for $ODE.

2.2.12 `$TABLE`

Syntax: C++
Multiple allowed: yes

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.17), 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 `$ERROR`

This is an alias for $TABLE.

2.2.14 `$EVENT`

Syntax: C++
Multiple allowed: no
mrgsolve version: >= 1.5.2

Use the $EVENT block to write code to implement doses or other modeled events inside your model file. $EVENT functions just like $ERROR (or $TABLE), but it gets called just prior to $ERROR. This lets you set up modeled doses or events for execution “now” and also seem them in the simulated output on that same record. For example

$EVENT 
if(NEWIND > 1) return;

evt::ev dose = evt::bolus(250, 2);
evt::ii(dose, 12);
evt::ss(dose, 1); 
evt::addl(dose, 9)

self.push(dose);

A subtle difference between $ERROR and $EVENT

Coding modeled events in $EVENT rather than $TABLE will give a very limited and subtle difference in output that you will most likely see when you administer a bolus dose to a compartment you are watching (e.g. a bolus to the central compartment). Using $EVENT will let you see the result of the dose immediately - on the record where the dose code is executed; when $TABLE is used, you will only see the result of that dose on subsequent output records. You’d see a difference between $TABLE and $EVENT behavior when you bolus into a depot compartment too, but you probably aren’t noticing because you are likely monitoring the central compartment, not the depot. Regardless, the $TABLE / $EVENT difference is just bookkeeping for the current record; you might see that bookkeeping difference on the current record, but the difference will not factor into how the system advances beyond that record.

The key point for using $EVENT

For the time being, the important thing to remember when you want to dose “now” (not some time in the future) is those interventions should be coded into $EVENT (or $ERROR) and they will happen right after the system advances in time on the current record, not prior to advancing in time for the next record. I recommend you always use $EVENT or $TABLE to execute doses or other events from inside your model. This distinction only applies to doses happening “now”; there is no difference in behavior when the dose is scheduled to happen at a future time.

2.2.15 `$PREAMBLE`

Syntax: C++
Multiple allowed: no

This is the fourth C++ code block. It is called once in two different settings (Chapter 8):

Immediately prior to starting the simulation run
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)

2.2.16 `$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 4.7 for additional information regarding data sets in use with $PRED block.

2.2.17 `$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.28).

Tip

In addition to listing model variables for output 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.18 `$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 6).

2.2.19 `$SIGMA`

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

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 EPS values.

The $SIGMA block functions like the $OMEGA block. See $OMEGA for details.

2.2.20 `$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.21 `$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:

Write #define preprocessor statements
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.

2.2.22 `$PKMODEL`

Syntax: R
Multiple allowed: no

This code block implements a one- or two-compartment PK model where the system is calculated by algebraic equations, not ODEs. mrgsolve handles the calculations and an error is generated if both $PKMODEL and $ODE blocks are included in the same model specification file.

This is an options-only block. The user must specify the number of compartments (1 or 2) to use in the model as well as whether or not to include a depot dosing compartment. See ?PKMODEL for more details about this block, including specific requirements for symbols that must be defined in the model specification file.

The $CMT or $INIT block must also be included with an appropriate number of compartments. Compartment names, however, may be determined by the user.

Example:

[ CMT ] GUT CENT PERIPH

[ PKMODEL ] ncmt=2, depot=TRUE

As of version 0.8.2, we can alternatively specify the compartments right in the $PKMODEL block:

$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.23 `$PLUGIN`

Syntax: text
Multiple allowed: yes

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

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.23.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.24 `$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.25 `$NMXML`

Syntax: R
Multiple allowed: yes

The $NMXML block lets you read and incorporate results from a NONMEM run into your mrgsolve model.

Syntax

If your model run is 1000 and the .xml file is located in 1000/1000.xml, use this syntax

$NMXML 
run = 1000
project = "../model/nonmem"
root = "cppfile"

If your model run is 1000 and the .xml file is located in some other location use the path argument to specify the complete path to the .xml file

$NMXML 
path = "../<path>/<to>/<output>/<files>/1000.xml"
root = "cppfile"

This latter syntax may be needed when using PsN and output files are organized in a way different from the previous example.

Details

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.18) and SIGMA will be captured as a $SIGMA block (Section 2.2.19). 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 path 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-455b1da05dc8 <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.26 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.26 `$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).

Syntax

If your model run is 1000 and the .ext file is located in 1000/1000.ext, use this syntax

$NMEXT 
run = 1000
project = "../model/nonmem"
root = "cppfile"

If your model run is 1000 and the .ext file is located in some other location use the path argument to specify the complete path to the .ext file

$NMEXT 
path = "../<path>/<to>/<output>/<files>/1000.ext"
root = "cppfile"

This latter syntax may be needed when using PsN and output files are organized in a way different from the previous example.

Details

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.

mod <- modlib("1005", compile = FALSE)

as.list(mod)$nm_import

See the $NMXML topic in Section 2.2.25 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.

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.27 `$INCLUDE`

Syntax: text
Multiple allowed: no

To include your own header file(s) in a model use $INCLUDE

$INCLUDE
mystuff.h
otherstuff.h

$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.28 `$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.23) 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.")  
}

$NMXML
run = 101
project = here::here("model/nonmem")

You can access the model environment using the env_get() function:

env_get(mod)

or directly in the envir slot

mod@envir

See also env_ls() to list the names of variables in the model environment, env_update() to change the values of objects in the model environment, or env_eval to re-evaluate the model environment code.

2.3 Variables and Macros

This section describes some macros and internal variables that can be used in model specification files. In the following section, we adopt the convention that CMT stands for a compartment in the model.

IMPORTANT NOTE:

It should be clear from the usage examples which variables can be set by the user and which are to be read or checked. All internal variables are pre-defined and pre-initialized by mrgsolve. The user should never try to declare an internal variable; this will always result in an compile-time error.

2.3.1 `ID`

The current subject identifier. ID is an alias for self.id.

2.3.2 `TIME`

Gives the time in the current data set record. This is usually only used in $MAIN or $TABLE. TIME is an alias for self.time. Contrast with SOLVERTIME.

2.3.3 `SOLVERTIME`

Gives the time of the current timestep taken by the solver. This is can only be used in $ODE. Contrast with TIME.

Starting with mrgsolve version 1.0.0, the variable T is provided as a synonym to SOLVERTIME inside $ODE when you invoke the nm-vars plugin (Section 10.2).

2.3.4 `T`

Only provided when the nm-vars plugin is invoked; see SOLVERTIME.

2.3.5 `EVID`

EVID is an event id indicator. mrgsolve recognized the following event IDs:

0 = an observation record
1 = a bolus or infusion dose
2 = other type event, with solver reset
3 = system reset
4 = system reset and dose
8 = replace

EVID is an alias for self.evid.

2.3.6 `CMT`

CMT is the current compartment number. In this case, CMT is used literally (not a stand-in for the name of a compartment). For example:

$CMT GUT CENT PERIPH

$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.7 `AMT`

AMT is the current value of dose amount.

AMT is an alias for self.amt.

2.3.8 `NEWIND`

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.9 `SS_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 9); 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.10 `simeta()`

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

$MAIN
simeta();

will re-simulate all ETA(n).

2.3.11 `simeps()`

simeps() works like simeta(), but all EPS(n) values are re-simulated rather than the ETA values.

2.3.12 `self`

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.13 `self.time`

The current data set TIME.

2.3.14 `self.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.15 `self.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.16 `self.nid`

The number of IDs in the data set.

2.3.17 `self.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.18 `self.nrow`

The total number of rows in the output data set.

2.3.19 `self.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.20 `self.envir`

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

For example

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

2.3.21 `self.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.22 `self.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.23 `self.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.

2.3.24 `self.stop()`

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.25 `self.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().

2.3.26 `self.stop_id_cf()`

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 current values are carried forward (cf) for the remaining output records for that individual.

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

2.3.27 `THETA(n)`

THETA(n) is ubiquitous in NONMEM control streams, representing estimated fixed effect parameters. mrgsolve does not attach any special meaning to THETA(n) but starting with version 1.0.0 it will translate THETA(n) to THETAn. This is useful when importing NONMEM estimates using $NMXML or $NMEXT, where THETAs are imported as parameters with names THETA1, THETA2, THETA3 and so on (see Section 2.2.25 and Section 2.2.26). THETA(n) is simply a NONMEM-oriented style for referring to THETAn.

2.3.28 `ETA(n)`

ETA(n) is the value of the subject-level variate drawn from the model OMEGA matrix. ETA(1) through ETA(25) have default values of zero so they may be used in a model even if appropriate OMEGA matrices have not been provided.

For example:

$OMEGA
1 2 3

$MAIN
double CL = TVCL*exp(ETA(1));
double VC = TVVC*exp(ETA(2));
double KA = TVKA*exp(ETA(3));

Here, we have a 3x3 OMEGA matrix. ETA(1), ETA(2), and ETA(3) will be populated with variates drawn from this matrix. ETA(4) through ETA(25) will be populated with zero.

2.3.29 `EPS(n)`

EPS(n) holds the current value of the observation-level random variates drawn from SIGMA. The basic setup is the same as detailed in ETA(n).

Example:

[ CMT ] CENT

[ PARAM ] CL=1, VC=20

[ SIGMA ] 
25 0.0025

[ TABLE ]
double DV = (CENT/VC)*(1+EPS(2)) + EPS(1);

2.3.30 `SIGMA(n)`

Starting with version 1.0.8, users have read-only access to on-diagonal elements of $\Sigma$ through the SIGMA() macro. For example

$SIGMA 0.1 12

$TABLE
double STD = sqrt(SIGMA(1) + pow(F,2)*SIGMA(2));

2.3.31 `table(<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.12 and also Section 2.2.17.

2.3.32 `F_CMT`

For the CMT compartment, sets the bioavailability fraction for that compartment.

Example:

$MAIN
F_CENT = 0.7;

2.3.33 `ALAG_CMT`

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

Example:

$MAIN
ALAG_GUT = 0.25;

2.3.34 `R_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.35 `Rn`

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

2.3.36 `D_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 10.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
until
now

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 CL = TVCL*exp(ETA(1));

We want CL 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 or $EVENT, 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 CL = TVCL*exp(ETA(1));

a double-precision numeric variable is created (CL) in the $MAIN block. When mrgsolve parses the model file, this code gets translated to

$GLOBAL
namespace {
  double CL;
}

$MAIN
CL = TVCL*exp(ETA(1));

That is, CL 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 CL variable in $TABLE:

$MAIN
double CL = TVCL*exp(ETA(1));
double V2 = TVV2*exp(ETA(2));

$TABLE
double KE = CL/V2;

$CAPTURE KE

To declare a variable that is local to a particular code block:

$MAIN

local_double CL = TVCL*exp(ETA(1));

The local_double 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.

Similarly, use local_int or local_bool to declare a variable that will not be made global by mrgsolve. One common application of local_int is when you need to write a for-loop

$TABLE 
for(local_int i = 0; i < 9; ++i) {}. 
  // do something with i
}

Here, we we want i declared locally for use in that for-loop and we use local_int to make that happen.

2.5.2 Initializing and resetting

mrgsolve takes a hands-off approach to variables created by the user; it will not initialize or reset variables from record to record or from individual to individual. Failing to properly initialize or reset user variables can have unintended consequences, especially when conditional logic is used.

Generally, a default value for all variables should be set prior to conditional assignment based on some test condition.

2.5.3 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;

[ 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.5.4 Automatic declaration

There is a plugin you can invoke which will instruct mrgsolve to scan your code and discover variables to be declared with type double. See the autodec plugin in Section 10.1 for details.

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 10 and Section 10.6).

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

2.1 How / where to write a model

2.1.1 Separate file

2.1.1.1 Using .cpp file extension

2.1.1.2 Using any file extension

2.1.1.3 The model project directory

2.1.2 Inline / code

2.1.3 Comments

2.2 Code blocks

2.2.1 About code blocks

2.2.2 Programmatic or bulk initialization

2.2.3 $PROB

2.2.4 $PARAM

2.2.5 $INPUT

2.2.6 $FIXED

2.2.7 $CMT and $INIT

2.2.8 $MAIN

2.2.9 $PK

2.2.10 $ODE

2.2.11 $DES

2.2.12 $TABLE

2.2.13 $ERROR

2.2.14 $EVENT

2.2.15 $PREAMBLE

2.2.16 $PRED

2.2.17 $CAPTURE

2.2.18 $OMEGA

2.2.19 $SIGMA

2.2.20 $SET

2.2.21 $GLOBAL

2.2.22 $PKMODEL

2.2.23 $PLUGIN

2.2.23.1 Some examples

2.2.24 $THETA

2.2.25 $NMXML

Syntax

2.2.26 $NMEXT

Syntax

Details

2.2.27 $INCLUDE

2.2.28 $ENV

2.3 Variables and Macros

2.3.1 ID

2.3.2 TIME

2.3.3 SOLVERTIME

2.3.4 T

2.3.5 EVID

2.3.6 CMT

2.3.7 AMT

2.3.8 NEWIND

2.3.9 SS_ADVANCE

2.3.10 simeta()

2.3.11 simeps()

2.3.12 self

2.3.13 self.time

2.3.14 self.cmt

2.3.15 self.amt

2.3.16 self.nid

2.3.17 self.idn

2.3.18 self.nrow

2.3.19 self.rown

2.3.20 self.envir

2.3.21 self.tad()

2.3.22 self.mtime(<time>)

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

2.3.24 self.stop()

2.3.25 self.stop_id()

2.3.26 self.stop_id_cf()

2.3.27 THETA(n)

2.3.28 ETA(n)

2.3.29 EPS(n)

2.3.30 SIGMA(n)

2.3.31 table(<name>)

2.3.32 F_CMT

2.3.33 ALAG_CMT

2.3.34 R_CMT

2.3.35 Rn

2.3.36 D_CMT

2.3.37 NONMEM-like syntax

2.4 Reserved words

2.5 Derive new variables

2.1.1.1 Using `.cpp` file extension

2.1.2 Inline / `code`

2.2.3 `$PROB`

2.2.4 `$PARAM`

2.2.5 `$INPUT`

2.2.6 `$FIXED`

2.2.7 `$CMT` and `$INIT`

2.2.8 `$MAIN`

2.2.9 `$PK`

2.2.10 `$ODE`

2.2.11 `$DES`

2.2.12 `$TABLE`

2.2.13 `$ERROR`

2.2.14 `$EVENT`

2.2.15 `$PREAMBLE`

2.2.16 `$PRED`

2.2.17 `$CAPTURE`

2.2.18 `$OMEGA`

2.2.19 `$SIGMA`

2.2.20 `$SET`

2.2.21 `$GLOBAL`

2.2.22 `$PKMODEL`

2.2.23 `$PLUGIN`

2.2.24 `$THETA`

2.2.25 `$NMXML`

2.2.26 `$NMEXT`

2.2.27 `$INCLUDE`

2.2.28 `$ENV`

2.3.1 `ID`

2.3.2 `TIME`

2.3.3 `SOLVERTIME`

2.3.4 `T`

2.3.5 `EVID`

2.3.6 `CMT`

2.3.7 `AMT`

2.3.8 `NEWIND`

2.3.9 `SS_ADVANCE`

2.3.10 `simeta()`

2.3.11 `simeps()`

2.3.12 `self`

2.3.13 `self.time`

2.3.14 `self.cmt`

2.3.15 `self.amt`

2.3.16 `self.nid`

2.3.17 `self.idn`

2.3.18 `self.nrow`

2.3.19 `self.rown`

2.3.20 `self.envir`

2.3.21 `self.tad()`

2.3.22 `self.mtime(<time>)`

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

2.3.24 `self.stop()`

2.3.25 `self.stop_id()`

2.3.26 `self.stop_id_cf()`

2.3.27 `THETA(n)`

2.3.28 `ETA(n)`

2.3.29 `EPS(n)`

2.3.30 `SIGMA(n)`

2.3.31 `table(<name>)`

2.3.32 `F_CMT`

2.3.33 `ALAG_CMT`

2.3.34 `R_CMT`

2.3.35 `Rn`

2.3.36 `D_CMT`

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