One Model, Several Software Representations
A pharmacometric model is more than a collection of software statements. It specifies how population-level quantities, individual differences, evolving system states, and measured data relate. NONMEM, Monolix, nlmixr2, and mrgsolve represent these roles using different project structures and syntax, but the underlying mathematical distinctions remain the same.
The purpose of a cross-tool map is not to make the four tools interchangeable. In particular, a tool used to estimate population parameters does not become an estimator simply because it can represent a parameter or simulate random variation. The useful question is whether each mathematical role has been represented consistently, and whether the tool is being used for a task it supports.
A compact mathematical description helps keep these roles separate. Let \(x_i(t)\) denote the state for individual \(i\), \(\theta\) the population-level structural parameters, and \(\eta_i\) that individual’s random effects. A general dynamic and observation model can be written as:
Here \(y_{ij}\) is the measured value for individual \(i\) at observation time \(t_{ij}\); \(\widehat y_{ij}\) is the model prediction; and \(\epsilon_{ij}\) is residual error. In a common hierarchical formulation, \(\eta_i\) has mean zero and covariance matrix \(\Omega\), while residual error has a specified distribution conditional on the model prediction. Normality and independence are modeling assumptions, not universal properties: for example, residuals may have prediction-dependent variance or a non-normal distribution. The chosen structure must match the data and estimation method.
The same symbol can appear under different software labels. A model’s interpretation comes from its equations and assumptions, not from the label alone. The broad correspondence is:
| Mathematical role | NONMEM | Monolix | nlmixr2 | mrgsolve |
|---|---|---|---|---|
| Population structural parameters | Typically population parameters such as THETA | Population parameters in the structural model | Population parameters in the model definition | Model parameters supplied for simulation |
| Dynamic state and evolution | Compartment or differential-equation representation | Structural model and its state equations | Model equations, including ODE-based components | Initial states and state-derivative equations |
| Predicted observation | Prediction and observation-model calculation, commonly involving F and Y | Structural prediction linked to the observation model | Model prediction linked to the residual-error specification | Calculated or captured model output |
| Interindividual random effects | ETA terms with an OMEGA covariance structure | Individual variability and its distribution | ETA terms and an associated covariance structure | Effects supplied or generated for simulation; not population estimation |
| Residual error | Observation model, often expressed in the Y calculation | Error model for observations | Residual-error model in the model definition | Must be represented in simulation output or added by the workflow |
Map the Equations Before the Syntax
When translating a model, first write down its roles independently of the software. Record the population parameters and units; the individual parameter relationships; the state equations and initial conditions; the function that produces a predicted observation; and the distribution of residual error. Then identify where each item is expressed in the target tool.
For example, a population parameter may be related to an individual parameter by an exponential random-effects model:
This says that \(k_{\mathrm{pop}}\) is the typical value on the parameter’s natural scale, \(\eta_i\) is a mean-zero random effect on the log scale, and \(\omega^2\) is the variance of that random effect. The individual parameter \(k_i\) is therefore positive. A software translation must preserve both the exponential relationship and the meaning of the variance. Supplying \(\omega\) where a tool expects \(\omega^2\), for instance, changes the model.
Worked Example: Translating a Population Parameter and an Individual Effect
Consider a hypothetical model with \(k_{\mathrm{pop}}=0.20\ \mathrm{h}^{-1}\), an individual random effect \(\eta_i=0.2231\), and the exponential relationship above. Calculate the individual parameter, then state what must remain the same when this relationship is represented in each tool.
Calculation. Substitute the given values:
The effect is dimensionless, so \(k_i\) retains the units \(\mathrm{h}^{-1}\). Across NONMEM, Monolix, and nlmixr2, the population parameter, random-effect distribution, and transformation must describe this same relationship. In mrgsolve, the same relationship can be used to calculate or simulate an individual parameter, but this does not estimate \(k_{\mathrm{pop}}\) or \(\omega^2\). Those estimates must come from an estimation workflow.
State variables must also be kept distinct from parameters. A parameter controls the evolution rule; a state is a quantity whose value changes over time according to that rule. A prediction is then calculated from the state and other model quantities. Some software may use a compartment-oriented representation, while another uses explicit state equations. Those are different representations of model structure, not evidence that a state variable and a parameter are interchangeable.
Worked Example: Separating a State, Prediction, and Residual
Suppose a hypothetical model produces a state value \(x_i(t)=8.0\ \mathrm{mg/L}\) at one observation time. For this illustration, the observation prediction is \(\widehat y_i=x_i(t)\), and the observation has additive residual error with standard deviation \(1.0\ \mathrm{mg/L}\). The measured value is \(y_i=9.2\ \mathrm{mg/L}\). Calculate the residual and its standardized value.
Residual. Using observed minus predicted:
Standardized residual. Dividing by the stated residual standard deviation gives \(1.2\ \mathrm{mg/L}/(1.0\ \mathrm{mg/L})=1.20\), a dimensionless value. The state is \(8.0\ \mathrm{mg/L}\); the prediction is \(8.0\ \mathrm{mg/L}\) because this example defines the observation function as the identity; and the residual is the difference from the measurement. A translation must preserve the observation function and additive error model, not merely calculate the same state.
In NONMEM, the prediction and observation model are represented in the model’s prediction/error calculation; Monolix separates structural prediction from the observation error model; nlmixr2 expresses prediction and residual error in its model specification. In mrgsolve, the model can calculate and return the prediction, while the simulation workflow must explicitly generate or add observation error if noisy observations are needed.
Random Effects, Residual Error, and Tool Capability
Random effects and residual error describe different sources of variation. Random effects represent modeled differences at a level such as individual; residual error represents observation-level deviations conditional on the model. Combining them into one unexplained noise term generally changes the model. An individual effect may persist across that person’s observations, whereas a separately generated residual can vary from observation to observation.
The distinction matters especially when interpreting mrgsolve alongside estimation tools. NONMEM, Monolix, and nlmixr2 are used for population-model estimation and can express population parameters, random-effects structures, and observation models within their estimation workflows. mrgsolve is principally a simulation tool. It can run a specified model with parameter values and can be used in workflows involving individual effects and simulated observation noise, but the presence of random-effect or error-related quantities in a simulation setup does not mean mrgsolve is estimating their population distributions from data.
Worked Example: Variance Is Not Standard Deviation
Suppose a hypothetical random effect satisfies \(\eta_i\sim N(0,0.09)\). An individual’s effect is \(\eta_i=0.15\), and \(k_{\mathrm{pop}}=0.20\ \mathrm{h}^{-1}\) in the exponential model. Find the random-effect standard deviation and the resulting individual parameter.
The second argument of \(N(\mu,\sigma^2)\) is the variance, so \(\omega^2=0.09\) and \(\omega=\sqrt{0.09}=0.300\). Substitution into the individual-parameter equation gives:
The value \(0.300\) is the standard deviation of \(\eta_i\), not its variance. A translation must check what quantity each input expects: a variance, standard deviation, or covariance matrix. Giving \(0.300\) as the variance would imply a different distribution, with standard deviation \(\sqrt{0.300}\approx0.5477\).
Why Equivalent Models Can Still Disagree
Two implementations are mathematically equivalent only when they encode the same model, not merely when their parameter names look similar. Numerical results may differ if units, initial conditions, transformations, random-effect distributions, or residual-error definitions differ. Even when these are aligned, numerical solvers and estimation settings can affect computed results. The comparison workflow should therefore establish model equivalence before attributing a discrepancy to estimation or optimization.
List the state equations, observation function, parameter relationships, and probability distributions, including units and scales.
Identify the corresponding parameter, state, prediction, random-effect, and residual-error representation in the target tool.
Use fixed parameter values and no random effects or residual noise to compare predictions at selected times. Then add individual effects and residual error separately.
Separate parameters estimated from data from values supplied for simulation, and record any differences in solver or estimation settings.
Useful diagnostics include checking prediction units and scales, plotting predicted values against time, and comparing the deterministic predictions before adding stochastic components. If deterministic predictions already differ, inspect state equations, initial conditions, parameter transformations, event handling, and units. If predictions match but noisy outputs differ, check random-number generation, random-effect covariance, residual-error variance, and whether error is generated once per individual or once per observation as intended.
A common mistake is to see an OMEGA-like declaration or a generated ETA in a simulation tool and conclude that it has fitted interindividual variability. That reasoning confuses a model component with an estimation procedure. The correction is to ask whether the tool is estimating the population distribution from observations or merely using supplied values or draws to generate simulated individuals. Likewise, a close match in one predicted curve is not proof that residual-error assumptions or random-effects covariance structures also match.
Limits of a Cross-Tool Map
The correspondence table is conceptual, not a guarantee of identical syntax, defaults, numerical behavior, or supported features. Each tool has its own conventions for data, model files, parameter declarations, estimation, and simulation. A translation can preserve the central equations while still requiring tool-specific choices for solver tolerances, event records, parameter constraints, and output handling. Those choices should be documented and tested rather than treated as invisible implementation details.
When comparing runs, keep inputs and model assumptions aligned before interpreting objective-function changes, parameter shifts, or precision differences. A discrepancy may reflect a changed model definition rather than a software-specific estimation outcome. Conversely, agreement in a small test case does not establish equivalence for all data, parameter values, or solver settings.
Check Your Understanding
Answer each question, then compare with the brief solutions below.
- A model predicts \(5.0\ \mathrm{mg/L}\), and the observed value is \(4.4\ \mathrm{mg/L}\). Under the observed-minus-predicted convention, what is the residual?
- If \(\eta_i\sim N(0,0.16)\), what is the standard deviation of the random effect?
- Why does representing individual effects in a mrgsolve simulation not by itself constitute estimation of a population random-effects distribution?
- Name two model components to check if deterministic predictions differ between tools before random effects or residual error are added.
Solutions. (1) \(4.4-5.0=-0.6\ \mathrm{mg/L}\). (2) The variance is \(0.16\), so the standard deviation is \(\sqrt{0.16}=0.400\). (3) Simulation uses supplied values or generated draws; population estimation requires fitting a distribution to observed data. (4) Examples include state equations, initial conditions, units, parameter transformations, and the observation function.