Tutorials › Pharmacometric Workflows and Reproducible Modeling › Transferring NONMEM Parameters into mrgsolve

Software Landscape and Interoperability · Tutorial 216 of 399

Transferring NONMEM Parameters into mrgsolve

You will transfer a NONMEM estimate set into mrgsolve while preserving the parameterization and distinguishing typical-subject predictions from population simulations.

Intermediate 11 min read

What You'll Learn

  • Identify which NONMEM estimates belong in mrgsolve fixed parameters and which describe random-effect variability.
  • Translate a NONMEM OMEGA variance-covariance matrix without confusing variances, standard deviations, and covariances.
  • Implement a lognormal clearance and volume model in mrgsolve.
  • Calculate a typical-subject prediction by setting individual random effects to zero.
  • Simulate a population using the transferred random-effect covariance structure.
  • Check the transfer for scale, ordering, units, and numerical consistency.

Why a Parameter Transfer Needs More Than a List of Estimates

A NONMEM run can provide estimates for typical parameters, interindividual variability, and residual error. Those estimates do not all enter a simulation in the same way. A typical clearance estimate is a fixed-effect value; an OMEGA entry describes the distribution of an individual random effect; and a SIGMA entry concerns residual variability. Treating all three as ordinary parameter values can produce a simulation that runs but represents a different model.

The aim here is to carry a fitted NONMEM parameter set into an mrgsolve model while retaining its meanings. This assumes the structural model and parameterization have already been mapped correctly, including units, covariate relationships, and event definitions. The earlier material on mapping parameter names and values is useful for that preliminary alignment; the focus here is on extracting estimates and using them for typical-subject and population simulation.

Definition: A fixed-effect estimate is an estimated population-level model quantity, such as a typical clearance. An OMEGA estimate is an element of the covariance matrix \(\Omega\) for interindividual random effects \(\eta\). A SIGMA estimate describes the residual-error random variable \(\epsilon\), under the residual model specified for the observation. These quantities have different roles and must be transferred according to their model definitions.

Read the Estimate in Its Model Context

Suppose a one-compartment IV bolus model uses lognormally distributed interindividual variability in clearance and volume:

$$ CL_i = TVCL\exp(\eta_{CL,i}), \qquad V_i = TVV\exp(\eta_{V,i}), \qquad \begin{pmatrix}\eta_{CL,i}\\ \eta_{V,i}\end{pmatrix} \sim N\left( \begin{pmatrix}0\\0\end{pmatrix}, \Omega \right). $$

Here \(TVCL\) and \(TVV\) are fixed effects, while \(\Omega\) is a variance-covariance matrix on the log-random-effect scale. With this parameterization, an individual with both random effects equal to zero has \(CL_i=TVCL\) and \(V_i=TVV\). For the lognormal parameters, those fixed effects are medians across the population, not generally arithmetic means.

A NONMEM report may show estimates as THETA values and OMEGA or SIGMA entries. Read the model code and the output labels together: THETA index alone does not tell you the parameter's units or whether a transformation is applied. Likewise, an OMEGA entry is not necessarily a standard deviation. In the model above, a diagonal OMEGA entry is a variance of an \(\eta\), and an off-diagonal entry is a covariance. If the output or a secondary export instead presents standard deviations or correlations, convert them back to the intended covariance matrix before entering that structure in the simulation model.

Worked Example: Extracting and Checking an Estimate Set

Consider a hypothetical NONMEM run for the model above. The final estimates are \(TVCL=4.00\ \mathrm{L/h}\), \(TVV=20.0\ \mathrm{L}\), and \[ \Omega = \begin{pmatrix} 0.0900 & 0.0120\\ 0.0120 & 0.0400 \end{pmatrix}. \] The rows and columns are ordered as \(\eta_{CL}\), then \(\eta_V\). Identify what belongs in the fixed parameters and check the covariance entries.

The fixed parameters are \(TVCL=4.00\ \mathrm{L/h}\) and \(TVV=20.0\ \mathrm{L}\). The diagonal entries are \(\operatorname{Var}(\eta_{CL})=0.0900\) and \(\operatorname{Var}(\eta_V)=0.0400\); the off-diagonal entries are \(\operatorname{Cov}(\eta_{CL},\eta_V)=0.0120\). The corresponding standard deviations are \(\sqrt{0.0900}=0.300\) and \(\sqrt{0.0400}=0.200\), but these standard deviations are not the values to enter when the model expects OMEGA variances.

The implied correlation is \[ \rho=\frac{0.0120}{\sqrt{0.0900}\sqrt{0.0400}} =\frac{0.0120}{0.300(0.200)}=0.200. \] The determinant is \(0.0900(0.0400)-(0.0120)^2=0.003456>0\), and both diagonal entries are positive. Thus this symmetric two-dimensional matrix is positive definite and is a valid covariance matrix for the stated random effects.

Map the Estimates into an mrgsolve Model

For this model, the fixed-effect values become named mrgsolve parameters, and the OMEGA estimates define the random-effect distribution. A compact model specification can be written as follows; the model code is illustrative and assumes the state is the amount in a central compartment.

Model specification: In the parameter block, set TVCL to 4.00 and TVV to 20.0, with units understood as L/h and L. Define an OMEGA block in the order \(\eta_{CL},\eta_V\), with rows \(0.0900\) and then \(0.0120,\ 0.0400\). In the main model block, define \(CL=TVCL\exp(ETA(1))\) and \(V=TVV\exp(ETA(2))\). The central amount follows \(dCENT/dt=-(CL/V)CENT\), and concentration is \(CENT/V\).

The order is essential: the first random effect in the model must correspond to the first row and column of \(\Omega\). If the OMEGA order is reversed without changing the parameter equations, the marginal variances may appear plausible while the covariance is attached to the wrong pair of parameters. Keep a written index map beside the model code, especially when the NONMEM run has multiple OMEGA blocks.

Only include a SIGMA estimate when the simulation is intended to include the observation-error process. For example, a proportional residual model has a different interpretation from an additive residual model, and its implementation must match the specified observation equation. Do not place a residual variance in OMEGA or use an OMEGA variance as if it were a residual variance. A deterministic concentration trajectory can be simulated without residual error when the target is the structural prediction rather than simulated observations.

Typical-Subject Simulation

A typical-subject prediction sets the individual random effects to their reference values, usually \(\eta_{CL}=\eta_V=0\). For lognormal parameter models, this makes the individual parameters equal to their typical values. It does not represent the average concentration over a heterogeneous population: nonlinear transformations mean the average of individual predictions need not equal the prediction at typical parameters.

For an IV bolus dose \(D\) into the central compartment, the amount after dose is \(D\exp[-(CL_i/V_i)t]\), and the concentration is that amount divided by \(V_i\). This follows by solving the first-order differential equation \(dCENT/dt=-(CL_i/V_i)CENT\) with initial amount \(CENT(0)=D\). The solution is \(CENT(t)=D\exp[-(CL_i/V_i)t]\), so \(C_i(t)=D/V_i\exp[-(CL_i/V_i)t]\).

Worked Example: Predicting a Typical Subject

Use the hypothetical estimates above, set both \(\eta\) values to zero, and give a \(100\ \mathrm{mg}\) IV bolus at time zero. Calculate the concentration at 8 h.

The typical parameters are \(CL=4.00\ \mathrm{L/h}\) and \(V=20.0\ \mathrm{L}\), so the elimination rate constant is \(CL/V=4.00/20.0=0.200\ \mathrm{h^{-1}}\). Substituting into the concentration equation gives \[ C(8)=\frac{100\ \mathrm{mg}}{20.0\ \mathrm{L}} \exp[-(0.200\ \mathrm{h^{-1}})(8\ \mathrm{h})] =5.00\ \mathrm{mg/L}\times\exp(-1.600) \approx 1.009\ \mathrm{mg/L}. \]

The result is the structural concentration for the typical subject under the stated bolus model. It does not include residual error, and it is not a population-average concentration.

Population Simulation and a Numerical Precision Check

To simulate a population, retain the fixed effects and draw each subject's random-effect vector from the transferred distribution, \(\eta_i\sim N(0,\Omega)\). In mrgsolve, this means the OMEGA block must contain the intended covariance matrix and the simulation must use it as the source of individual random effects. A typical-subject run instead uses zero random effects; do not assume that supplying OMEGA alone makes every subject typical. Check the simulation settings and whether random effects are drawn or supplied as data.

The effect of a draw is applied through the parameter equations. For example, if a sampled subject has \(\eta_{CL}=-0.200\) and \(\eta_V=0.100\), then \(CL_i=4.00\exp(-0.200)=3.275\ \mathrm{L/h}\), and \(V_i=20.0\exp(0.100)=22.10\ \mathrm{L}\). Its concentration after a \(100\ \mathrm{mg}\) bolus at 8 h is \[ C_i(8)=\frac{100}{22.10}\exp\left[-\left(\frac{3.275}{22.10}\right)8\right] \approx 1.383\ \mathrm{mg/L}. \] This is one hypothetical subject conditional on the stated random-effect values, not an estimate of a population summary.

Calculation precision matters when reporting simulated concentrations. For instance, if the unrounded parameter formulas give \(D/V_i\approx2.7634\ \mathrm{mg/L}\) and \(\exp[-(CL_i/V_i)8]\approx0.25925\), then their product is \(2.7634(0.25925)=0.716411\ldots\), which rounds to \(C_i(8)\approx0.7164\ \mathrm{mg/L}\). If the factors are instead rounded to \(2.763\) and \(0.2592\) before multiplication, their product is \(0.7161696\), or \(0.7162\ \mathrm{mg/L}\). Use full precision for the calculation and round only the final reported result.

Derivation and Checks for the Random-Effect Transfer

A covariance matrix must be symmetric and positive semidefinite: for every vector \(a\), \(a^\mathsf{T}\Omega a\geq0\). This condition ensures that \(\Omega\) can describe a multivariate distribution and can be used to generate random effects. For a positive-definite matrix, a Cholesky factorization gives \(\Omega=LL^\mathsf{T}\). If \(z\) is a vector of independent standard normal draws, then \(\eta=Lz\) has mean zero and covariance \[ \operatorname{Cov}(\eta)=\operatorname{Cov}(Lz) =L\operatorname{Cov}(z)L^\mathsf{T} =LI L^\mathsf{T} =\Omega. \] This is why preserving the complete matrix, including its off-diagonal entries and ordering, matters for population simulation.

For the lognormal model, \(\eta_{CL}\sim N(0,\omega^2_{CL})\) implies \(\log(CL_i)=\log(TVCL)+\eta_{CL}\). Its median is \(TVCL\), because the median of a zero-mean normal random variable is zero. Its arithmetic mean is \(TVCL\exp(\omega^2_{CL}/2)\). With \(TVCL=4.00\ \mathrm{L/h}\) and \(\omega^2_{CL}=0.0900\), the arithmetic mean is \(4.00\exp(0.0450)\approx4.184\ \mathrm{L/h}\). Thus the typical value and the arithmetic mean are distinct, even though both can be called population summaries in different contexts.

Diagnostics, Alternatives, and Common Failure Modes

Before relying on a transferred simulation, inspect the inputs and outputs. Confirm that the final estimates—not initial values or an intermediate run—were copied; that the THETA mapping reflects the actual model equations; that units agree; and that OMEGA order, diagonals, and covariances are preserved. Check that the matrix is symmetric and positive semidefinite. Also verify the random-effect source: a simulation that inadvertently uses zero \(\eta\) values will reproduce a typical subject even if OMEGA is present.

A useful numerical diagnostic is to compare a small hand calculation with the simulation for a known subject. Set the random effects to zero, administer a single bolus, and compare the concentration at a selected time with the analytical expression. Then test a supplied nonzero random-effect vector and verify that the resulting clearance and volume match their exponential equations. These checks reveal many scale, index, and event-setup errors before a large population run.

The main alternatives are to transfer the estimated OMEGA covariance matrix and let mrgsolve generate random effects, or to generate and supply \(\eta\) values externally. The first is concise and keeps the distribution in the model; the second can be useful for controlled draws, reproducible comparisons, or nonstandard sampling schemes. In either case, the random effects must follow the intended distribution and use the same ordering and transformation. Neither alternative repairs a mismatch in model parameterization.

  • Confusing variance with standard deviation. If the model expects an OMEGA variance of \(0.0900\), entering \(0.300\) substitutes its standard deviation and changes the variability.
  • Dropping covariance. Replacing the example's \(0.0120\) covariance with zero changes the joint distribution, even though each marginal variance stays the same.
  • Calling the zero-random-effect prediction the population mean. It is the typical-subject prediction under this parameterization; nonlinearity generally makes it different from the average over simulated subjects.
  • Rounding intermediate factors. Premature rounding can alter the last reported digit, as the concentration calculation above shows.

These transfers are only as reliable as the correspondence between the fitted and simulation models. If the structural equations, transformations, covariates, or residual model differ, copied estimates alone cannot guarantee that the simulation represents the fitted model. The task here is faithful parameter transfer, not a substitute for checking the model definitions or the assumptions behind the estimation.

Check Your Understanding

Answer each question, then compare with the solutions below.

  1. For lognormal clearance with \(TVCL=5.00\ \mathrm{L/h}\) and \(\operatorname{Var}(\eta_{CL})=0.160\), what is the clearance at \(\eta_{CL}=0\), and what is the population arithmetic mean under the model?
  2. A report gives an OMEGA standard deviation of \(0.250\) for one random effect. What variance belongs in a model specification that expects OMEGA variance?
  3. Why does reversing the order of two random effects without changing their parameter equations cause a transfer error when OMEGA has a nonzero covariance?
  4. For the example's matrix, calculate the correlation between clearance and volume random effects.

Solutions. (1) At \(\eta=0\), clearance is \(5.00\ \mathrm{L/h}\). The arithmetic mean is \(5.00\exp(0.160/2)=5.00\exp(0.0800)\approx5.416\ \mathrm{L/h}\). (2) The variance is \(0.250^2=0.0625\). (3) The covariance is associated with matrix row and column order; reversing indices assigns the joint dependence to the wrong random-effect pair unless the equations and matrix are permuted together. (4) The correlation is \(0.0120/(\sqrt{0.0900}\sqrt{0.0400})=0.0120/(0.300\times0.200)=0.200\), dimensionless.

Clinical Biostats