From the Pharmacokinetic Model to a NONMEM Control Stream
A NONMEM control stream is an executable description of a model, the data it uses, and the estimation and reporting tasks to perform. For a one-compartment intravenous bolus model, the key implementation choices are how clearance and volume determine elimination, how the administered dose enters the compartment, how model amounts become concentrations, and how between-individual and residual variability are represented. A control stream can be syntactically valid while expressing the wrong model, so each piece should be checked against its mathematical meaning.
Assume a dose \(D\) is given as an instantaneous intravenous bolus at time zero. The body is represented by one well-mixed compartment with volume \(V\), and drug leaves that compartment by first-order elimination with clearance \(CL\). The model assumes \(CL>0\), \(V>0\), and constant values of both parameters over the modeled period. With amount \(A(t)\), the elimination rate constant is \(k=CL/V\). The concentration is \(C(t)=A(t)/V\).
This deterministic equation describes the typical structural model. A population analysis can allow individual \(CL\) and \(V\) to vary around typical values and can also allow measured concentrations to differ from their individual predictions. The control stream needs to keep those sources of variation distinct.
A Minimal Control Stream and What Each Part Does
The following hypothetical stream uses a comma-delimited data file with a header row. It expects one row per dose or observation event. The details of the input columns must match the file actually supplied to NONMEM.
Worked Example: Building the Core NONMEM Stream
Suppose a hypothetical dataset contains subject identifier, time in hours, observed concentration in milligrams per liter, dose amount in milligrams, event identifier, and a missing-dependent-variable indicator. A compact model specification is:
The section names begin with a dollar sign in a NONMEM control stream. Here they are shown in display math so the full stream is visually grouped. $INPUT maps the file columns to names used by the model; $DATA identifies the file and tells NONMEM to ignore its header row. $SUBROUTINES selects the built-in one-compartment bolus model and its clearance-volume parameterization.
In $PK, the THETA values are typical clearance and volume: the initial values are 5 L/h and 50 L. The exponential ETA terms keep individual parameter values positive. The two $OMEGA entries initialize the variances of the corresponding ETA terms: 0.09 for clearance and 0.04 for volume. These are variances on the log scale, not variances in L/h or L. S1=V scales the compartment amount to concentration for the dependent-variable prediction. In $ERROR, the proportional residual model adds an error term whose size is proportional to the prediction. $SIGMA initializes the variance of that dimensionless residual error to 0.04.
The estimation line requests first-order conditional estimation with the interaction option. The covariance step requests an estimate of parameter uncertainty when the fit supports it; it is not a guarantee that an uncertainty estimate will be available or reliable. The table requests selected data and model quantities for inspection. Exact table contents can depend on the NONMEM version and the variables made available by the model.
The stream is only consistent if the file matches its declarations. For example, a dose record should have a valid ID, TIME, and AMT, with EVID identifying it as a dose event and MDV indicating that there is no observed DV on that row. An observation record should have an observed DV, usually AMT equal to zero, and event and missingness indicators appropriate to an observation. The row conventions and values must be checked in the actual input file; the control stream cannot repair a misclassified record.
How the Structural Parameters Determine Predictions
The ADVAN1/TRANS2 choice lets the model be written in terms of clearance and volume rather than requiring the user to code the amount differential equation. The relationship \(k=CL/V\) is still essential for checking whether the parameter values and predicted concentrations make sense. The model solution follows by separating variables and integrating from the post-dose amount \(D\).
For \(A(t)>0\), \(dA/A=-(CL/V)dt\).
\(\ln[A(t)/D]=-(CL/V)t\), using \(A(0^+)=D\).
\(A(t)=D\exp[-(CL/V)t]\), so \(C(t)=D/V\exp[-(CL/V)t]\).
This derivation assumes a single instantaneous dose, a constant volume, and first-order elimination. It does not describe an infusion, absorption after an extravascular dose, time-varying clearance, or multiple compartments. Those situations require a different model specification or event structure.
Worked Example: Checking a Predicted Concentration
For a hypothetical individual, let \(D=500\) mg, \(CL=5\) L/h, and \(V=50\) L. Calculate the elimination rate constant and the concentration two hours after the bolus.
First, \(k=CL/V=(5\ \mathrm{L/h})/(50\ \mathrm{L})=0.1000\ \mathrm{h}^{-1}\). The post-dose concentration is \(D/V=(500\ \mathrm{mg})/(50\ \mathrm{L})=10.00\ \mathrm{mg/L}\). Therefore:
The concentration decreases from 10.00 mg/L immediately after the bolus to about 8.187 mg/L at two hours. A model output far from this value, given these same individual parameters and event time, would prompt checks of dose units, time units, parameter values, and concentration scaling.
Interindividual Variability and Residual Error
A common population parameterization assigns each individual a typical value multiplied by an exponential random effect. For example, \(CL_i=TVCL\exp(\eta_{CL,i})\) and \(V_i=TVV\exp(\eta_{V,i})\), where \(TVCL\) and \(TVV\) are typical values. The standard assumption for this simple specification is that the ETA terms have mean zero and are normally distributed; a diagonal OMEGA specification additionally assumes zero covariance between the two ETA terms. The exponential form ensures positive individual parameters. With mean-zero normal ETA, the typical value is the median of the corresponding lognormal individual-parameter distribution.
The residual model in the example is \(Y=IPRED(1+EPS)\), where \(IPRED\) is the individual prediction and \(EPS\) is a mean-zero normal residual term with variance SIGMA. Conditional on the individual prediction, this model has expected observation \(IPRED\) and residual standard deviation \(|IPRED|\sqrt{SIGMA}\). The proportional model is most directly suitable when variability grows roughly with concentration. It can be problematic near zero, where its implied residual standard deviation also approaches zero.
Worked Example: Interpreting ETA and SIGMA Values
Use the hypothetical typical parameters \(TVCL=5\) L/h and \(TVV=50\) L. Suppose an individual has \(\eta_{CL}=0.20\) and \(\eta_V=-0.10\). Calculate the individual parameters, then calculate the residual standard deviation at an individual prediction of 8.00 mg/L when SIGMA is 0.04.
The individual clearance is \(CL_i=5\exp(0.20)=6.107\) L/h, and the individual volume is \(V_i=50\exp(-0.10)=45.24\) L. Both remain positive, as required by the exponential parameterization. For the residual model, \(\sqrt{SIGMA}=\sqrt{0.04}=0.2000\), so the residual standard deviation at 8.00 mg/L is \(8.00(0.2000)=1.600\) mg/L. The conditional observation variance is therefore \(1.600^2=2.560\ (\mathrm{mg/L})^2\).
This calculation interprets SIGMA as a residual variance, not as a residual standard deviation. It also shows why the proportional-error coefficient is dimensionless: the predicted concentration supplies the units of the residual standard deviation.
Diagnostics, Alternatives, and Limits
Before interpreting estimates, check that the implemented model matches the intended design and units. A small set of targeted checks often catches problems that an estimation run alone cannot identify.
- Check event records. Confirm that dose rows carry the intended amount and event type, and that observation rows carry DV and are not accidentally marked missing. Verify that time is consistently in hours if CL is in L/h.
- Check concentration scaling. With S1=V and DV in mg/L, the model converts amount in mg to concentration in mg/L. Omitting or mis-specifying the scale can make predictions incompatible with the observed DV units.
- Check the implied curve. For a single bolus, concentrations should decline exponentially after the dose. The initial concentration is \(D/V\), and the rate of decline depends on \(CL/V\). Compare individual predictions with these expectations.
- Check parameter estimates and uncertainty. Implausible values, very large standard errors, or an unsuccessful covariance step can indicate weak information, poor starting values, an overly complex variability structure, or a data/model mismatch. They do not identify a unique cause by themselves.
- Check residual patterns. Trends in residuals over time or concentration suggest that the structural or residual-error assumptions may not describe the data adequately. A proportional model may be unsuitable when measurement error has a substantial concentration-independent component.
An alternative is to code the amount equation explicitly in a differential-equation block, such as a general ADVAN with user-specified differential equations. That can be useful when the structural model is not available as a built-in analytic model, but it adds implementation choices that must be checked, including state definitions and solver settings. For this simple bolus model, ADVAN1/TRANS2 expresses the intended structure more compactly. Another alternative is an additive or combined residual-error model; those address different patterns of measurement variability and should be selected based on the data and scientific assumptions, not merely because a model fails to converge.
The one-compartment model is a simplification: it assumes immediate distribution into one homogeneous compartment and a single first-order elimination process. A good fit does not establish that these assumptions are biologically exact. Nor does a successful estimation run prove the data were parsed or parameterized as intended. Reproducible modeling therefore means preserving the control stream, input data, run settings, and output together, and recording any changes between runs.
Common Mistakes
Correction: Units must be connected explicitly. In this model, setting \(S1=V\) provides the amount-to-concentration scaling, provided the volume and DV units are consistent. Check a hand-calculated prediction, such as \(D/V\) immediately after the bolus, against the model output.
Other frequent errors include treating an OMEGA entry as an ETA standard deviation rather than a variance, using a residual variance as though it were a residual standard deviation, or mixing hours and minutes between TIME and CL. Each can change the interpretation or numerical behavior of the model even if the control stream runs.
Check Your Understanding
Answer the questions, then compare your reasoning with the solutions.
- For \(CL=4\) L/h and \(V=40\) L, what is \(k\), and what is the half-life under this model?
- For a 400 mg bolus with \(V=40\) L, what is the concentration immediately after the dose?
- If the clearance ETA variance is 0.09, what is its ETA standard deviation? Is 0.09 measured in L/h?
- At an individual prediction of 5.00 mg/L and SIGMA=0.04, what residual standard deviation does the proportional-error model imply?
- Why should a modeler inspect dose and observation records even after NONMEM completes an estimation run?
Solutions. (1) \(k=4/40=0.1000\ \mathrm{h}^{-1}\). The half-life is \(\ln(2)/k=0.6931/0.1000=6.931\) h. (2) \(C(0^+)=D/V=400/40=10.00\) mg/L. (3) The ETA standard deviation is \(\sqrt{0.09}=0.3000\); the variance is on the log scale and is not measured in L/h. (4) The residual standard deviation is \(5.00\sqrt{0.04}=1.000\) mg/L. (5) A completed run does not establish that records were assigned the intended event, dose, observation, and missingness meanings; incorrect records can yield estimates for the wrong analysis.