mrgsolve Model Files and R Interfaces

From Model Description to a Runnable mrgsolve File

A pharmacometric model is not only a set of equations. It must also specify what quantities the software treats as parameters, how its dynamic states begin, how those states change over time, and which calculated quantities should be returned to the analyst. In mrgsolve, these tasks are organized into named blocks in a model file, which is compiled and controlled through R.

This tutorial focuses on the file structure and the R interface. The example is deliberately small: a single central compartment after an intravenous bolus. Its purpose is to show how model components connect and how to check a run, not to develop an estimation workflow. The model assumes a fixed dose, fixed positive clearance and volume, first-order elimination, and no absorption, random effects, or residual error.

The block names are part of the model’s interface, not decorative headings. For this example, $PARAM supplies clearance and volume, $INIT sets the starting central amount, $ODE defines its time derivative, and $TABLE with $CAPTURE defines and returns plasma concentration. A $PROB block can label the model for a human reader. A compact file can omit blocks that are not needed.

States, Parameters, and the Differential Equation

Let \(A(t)\) denote the amount of drug in the central compartment at time \(t\), measured in mg. Let \(CL\) be clearance in L/h and \(V\) be volume in L. With linear elimination, the rate of loss of amount is clearance multiplied by concentration. Since concentration is \(A(t)/V\), the differential equation is

The ratio \(CL/V\) has units \(h^{-1}\), so multiplying it by an amount in mg gives a rate in mg/h, as required for the derivative of amount with respect to time. The concentration reported at any time is \(CP(t)=A(t)/V\), in mg/L. These units are part of the model specification: if a dose is entered in mg but volume is supplied in mL without conversion, the resulting concentration will be wrong by a factor of 1,000.

A minimal model file can be written as the following sequence of blocks. The block contents are shown as separate lines; they belong together in one model source file.

$PROB One-compartment bolus demonstration

$PARAM CL = 5, V = 50

$INIT CENT = 0

$ODE

dxdt_CENT = -(CL/V)*CENT;

$TABLE

double CP = CENT/V;

$CAPTURE CP

Here CENT is the state corresponding to \(A(t)\). The derivative name begins with dxdt_ followed by the state name, so dxdt_CENT specifies the rate of change of CENT. The parameter values 5 L/h and 50 L imply an elimination rate constant of \(5/50=0.1\ h^{-1}\). The initial condition is zero because the system is empty before the dose. The dose is supplied as an event from R, rather than written into the initial-condition block.

The $PARAM values are defaults that can be changed when configuring a model in R. They are not automatically estimates, covariates, or random effects. The $TABLE block calculates a derived quantity after the state has been updated for the relevant time, and $CAPTURE makes CP available as a result column. Capturing a quantity does not alter the differential equation.

Compiling and Simulating From R

Save the model text in a source file named pk1.cpp in the working project directory. With mrgsolve installed, load the package and read the model by its name without the file extension. The following R commands compile the model and request output every hour through 24 hours:

library(mrgsolve)

mod <- mread("pk1", project = ".")

out <- mrgsim(mod, events = ev(amt = 100, cmt = 1), end = 24, delta = 1)

result <- as.data.frame(out)

In this event, amt is the dose amount in the model’s amount units and cmt = 1 identifies the first declared state, CENT. The event is an intravenous bolus. The simulation horizon is 24 hours, and delta = 1 requests hourly output. The resulting data frame includes time, the state CENT, and the captured quantity CP. A dose event changes the state; it is not a parameter declaration or a new differential equation.

The model can also be configured with new parameter values in R. For example, param(mod, CL = 6) creates a model configuration with clearance changed to 6 while retaining the other parameter value. The event and simulation settings can then be supplied to mrgsim again. Keep a record of which model configuration and event data generated each result; otherwise, a plausible output can still be associated with the wrong settings.

Numerical Checks and Failure Modes

A simulation is useful only if its output is consistent with the source model and the intended inputs. For this one-compartment bolus example, several checks follow directly from the stated assumptions. After the dose, CENT should be nonnegative and decline monotonically because \(dCENT/dt=-(CL/V)CENT\) is negative whenever \(CL>0\), \(V>0\), and CENT is positive. CP should equal CENT divided by V at every output time. If clearance is zero, the state should remain constant after the bolus; if volume increases with clearance unchanged, the initial concentration should decrease.

Common failures have recognizable patterns. A compiler error often means a typo, missing declaration, or malformed block. A run that completes but has no meaningful dose response can result from an event assigned to the wrong compartment. Concentrations off by a constant factor often indicate a unit mismatch or a mistaken volume. A flat profile despite positive clearance can signal that the state derivative was omitted or named incorrectly. Output lacking CP may mean the calculated variable was not captured. Solver warnings or unstable values call for checking the equations and units, parameter ranges, and integration settings rather than simply accepting the output.

Choosing the Right Level of Model Complexity

A useful alternative to a source model is a model assembled through an R interface or a higher-level modeling language. Such tools can be more convenient for defining estimation problems, data mappings, and model structures. mrgsolve is particularly direct when the goal is to specify a dynamic system and simulate it with explicit model code. The best choice depends on the task; the syntax is not itself evidence that one model is more scientifically appropriate than another.

The compact example also makes its limitations clear. It represents a single bolus and a single amount state, with fixed parameters and no measurement error. It does not establish how a dataset should be prepared, how parameters should be estimated, or how more complex dosing records and model features should be translated. Those additions change the model and its assumptions, so they should be introduced only when required and checked individually.