Model Compilation and Execution in mrgsolve

From Model Specification to Simulation

A model written for mrgsolve is a specification of state variables, parameters, differential equations, and calculated outputs. Before it can simulate a regimen, mrgsolve turns the model code into compiled code that can be called from R. The simulation then supplies values and events to that executable model and numerically advances the system over time.

This distinction between model code and run inputs is operationally important. Changing an ODE changes what the model does and requires the model code to be compiled again. Changing a clearance value or dose, while leaving the model structure unchanged, is ordinarily a runtime change: the same compiled model can be run again with new inputs.

Definition: Compilation is the translation of a model specification into executable code. Execution is a simulation run using that executable code, together with parameter values, event records, and requested output times. A compiled model can be executed repeatedly with different run inputs, provided those inputs are compatible with the compiled model.

In an mrgsolve model, sections such as $PARAM, $CMT, $ODE, $TABLE, and $CAPTURE describe parts of the model. The compiler uses this specification to produce code that mrgsolve can execute. R is used to create or read the model, compile it, set or update inputs, request simulation, and inspect the resulting output. Compilation is not parameter estimation: it does not find parameter values that best fit data.

A Small Model and Its Compilation

Consider a hypothetical one-compartment model for an intravenous bolus. Let $A(t)$ be the amount of drug in the compartment, $CL$ clearance, and $V$ volume of distribution. The model assumes that elimination is first order, with the rate of elimination equal to $CL$ times concentration. Since concentration is $A(t)/V$, the amount changes according to the following equation.

\frac{dA(t)}{dt}=-CL\frac{A(t)}{V}

The model treats $CL$ and $V$ as constant during a run, uses one well-mixed compartment, and represents the bolus as an instantaneous increase in amount. These are simplifying assumptions for this example, not claims about any particular medicine.

Worked Example: Compile and Run a One-Compartment Model

The following R code defines the model, compiles it with mcode(), and requests concentrations every 6 hours through 24 hours. The hypothetical parameters are $CL=5$ L/h and $V=50$ L, and the intravenous bolus is 100 mg.

R model specification and run: model_code <- '$PARAM CL=5, V=50\n$CMT CENT\n$ODE\ndxdt_CENT = -(CL/V)*CENT;\n$TABLE\ndouble CP = CENT/V;\n$CAPTURE CP'; mod <- mcode("onecomp", model_code); out <- mrgsim(mod, events=ev(amt=100, cmt=1), tgrid=seq(0, 24, by=6))

Here CENT is the amount state, and cmt=1 sends the dose to the first (and only) compartment. CP is calculated as amount divided by volume and captured in the output. The parameter values in $PARAM are defaults for the run; they are not inferred from the dose or from observations.

The model has elimination rate constant $k=CL/V=5/50=0.1000$ h$^{-1}$. Immediately after the bolus, $A(0)=100$ mg, so concentration is $100/50=2.000$ mg/L. The analytic concentration at time $t$ is $2.000e^{-0.1000t}$ mg/L. For example, at 6 hours it is $2.000e^{-0.6000}=1.098$ mg/L, and at 12 hours it is $2.000e^{-1.200}=0.6024$ mg/L. The simulated values should be close to these predictions, subject to numerical integration tolerance and output conventions.

Why the Analytic Solution Is a Useful Check

The differential equation also provides a simple independent check on execution. With constant $CL$ and $V$, write $k=CL/V$; then $dA/dt=-kA$. Separating variables and integrating from time zero to time $t$ gives the amount solution directly.

\frac{dA}{A}=-k\,dt,\qquad \int_{A(0)}^{A(t)}\frac{dA}{A}=-k\int_0^t dt,\qquad \log\!\left(\frac{A(t)}{A(0)}\right)=-kt

A(t)=A(0)e^{-kt},\qquad C(t)=\frac{A(t)}{V}=\frac{\text{Dose}}{V}e^{-(CL/V)t}

The derivation assumes a single compartment, constant positive $CL$ and $V$, and no additional input after the bolus. It does not describe oral absorption, time-varying clearance, or multiple compartments. Under these stated conditions, comparing simulated values with the analytic curve can reveal a wrong parameter, wrong compartment assignment, unit error, or unexpected dose timing. Small differences from numerical integration are expected; large or systematic differences call for investigation.

Runtime Parameter Updates and Repeated Runs

A parameter value can be changed on the model object without changing the equations. In mrgsolve, param() updates parameter values for subsequent execution. This is different from editing the model source: an updated clearance changes the rate used by the existing equation, while a change to the equation itself changes the compiled model specification.

Worked Example: Change Clearance Without Changing the Model

Use the same compiled one-compartment model and the same 100 mg bolus, but change clearance from 5 to 10 L/h while keeping volume at 50 L. A new run can be requested as follows: mod_hiCL <- param(mod, CL=10); out_hiCL <- mrgsim(mod_hiCL, events=ev(amt=100, cmt=1), tgrid=c(0, 12)).

The new rate constant is $k=10/50=0.2000$ h$^{-1}$. At 12 hours, the analytic prediction is $C(12)=2.000e^{-(0.2000)(12)}=2.000e^{-2.4}=0.1814$ mg/L. With the original clearance, the prediction was $2.000e^{-1.2}=0.6024$ mg/L. The lower concentration in the updated run follows from the higher clearance under this model; it is not a result of recompiling a different equation.

For repeated scenario runs, keep the compiled model and vary only the intended inputs. For example, the following code requests runs with two hypothetical bolus amounts: doses <- c(100, 200); dose_runs <- lapply(doses, function(d) mrgsim(mod, events=ev(amt=d, cmt=1), tgrid=c(0, 12))). At 12 hours, the model predicts 0.6024 mg/L for 100 mg and 1.205 mg/L for 200 mg. The second value doubles because the model is linear in dose when $CL$ and $V$ are fixed.

Each call to mrgsim() requests a simulation; it does not mean that the previous run's final amount is automatically carried into the next run. The event schedule and initial conditions determine the run. When simulating multiple individuals, individual-level parameter inputs can likewise be supplied to a run, while the model equations remain shared. Keep track of which values are defaults, which are overridden for a run, and which differ by individual.

A Practical Execution Workflow

A reproducible workflow separates model definition, compilation, and run specification. This makes it easier to establish whether an unexpected result came from the equations, the chosen inputs, or the execution request.

Specify the model.
Define parameters, compartments, equations, and captured outputs. Check that the dimensions are compatible; for example, $CL/V$ has units of inverse time.

Compile the specification.
Use an mrgsolve compilation workflow such as mcode() for model code supplied in R, or load and compile a model from a file. Compilation requires an available, working C++ toolchain.

Set run inputs.
Supply parameter values, event records, and requested output times. Confirm units, dose compartment, dose time, and initial conditions before execution.

Execute and check.
Call the simulation, inspect the expected output columns and time points, and compare a simple case with an analytic result or another justified reference.

Record what changed.
Retain the model source, run inputs, software and toolchain versions, and relevant settings so that a later run can be reconstructed and compared.

Diagnostics, Alternatives, and Limits

Compilation errors and simulation problems are different failure modes. A compilation error generally means the model source or compilation environment could not produce executable code; inspect the compiler message and the relevant model section. A successful compilation does not establish that the model is scientifically appropriate or that its units and event inputs are correct. Execution can still produce unexpected results because a dose was assigned to the wrong compartment, a parameter was not updated as intended, or output was requested at the wrong times.

Several checks are especially useful:

Confirm that the compiled model corresponds to the source version you intended to run, particularly after editing model code.
Print or inspect parameter values after an update, and verify that the values have the units expected by the equations.
Check the event time, amount, and compartment, then confirm that the requested time grid includes the times you need.
Look for missing, non-finite, or implausible output values and compare a simple case against a hand calculation where possible.
If a small code change alters model structure, compile the updated specification rather than assuming a runtime parameter update can change the equation.

Compiling model code and then changing runtime inputs is usually more efficient and clearer than repeatedly rebuilding the model for scenarios that differ only in parameters or events. Conversely, editing an ODE or adding a state or calculated variable changes the specification and calls for a new compilation step. A separate alternative is to calculate a simple model analytically, as above; this is valuable as a check but cannot replace execution for complex systems. Numerical simulations also have limitations: results depend on the model assumptions, the supplied inputs, and numerical solver settings. A compiled model is executable code, not evidence that those assumptions are adequate.

Common pitfall: Treating every change as either “just a parameter update” or “a new model” without distinguishing model structure from run inputs. Changing $CL$ from 5 to 10 L/h changes a runtime value in the example; changing the ODE to include another elimination pathway changes the model specification. The former can use the existing compiled model, while the latter requires updating and compiling the changed code.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions.

For the example model, what are the rate constant and predicted concentration at 6 hours when $CL=5$ L/h, $V=50$ L, and dose is 100 mg?
Does changing $V$ from 50 L to 60 L alter the differential equation's structure or a runtime parameter value? What should be checked before the next run?
Why can a model compile successfully and still produce an incorrect concentration-time profile?
Would changing the right-hand side of the ODE from $-CL \times CENT/V$ to a model with an additional elimination pathway be handled by updating $PARAM alone?

Solutions. (1) $k=5/50=0.1000$ h$^{-1}$, and $C(6)=2.000e^{-0.6}=1.098$ mg/L. (2) It changes a runtime parameter value, not the equation structure. Check that the updated value is used and that the units remain consistent; the new initial concentration after a 100 mg bolus would be $100/60=1.667$ mg/L. (3) Compilation checks whether executable code can be produced, not whether the dose, parameter units, event timing, compartment assignment, or scientific assumptions are correct. (4) No. An additional pathway changes the model equations, so the model source must be changed and the updated specification compiled.