A One-Compartment Model in nlmixr2

From the Pharmacokinetic Model to nlmixr2

A software model is a precise statement of a pharmacokinetic model, not merely a collection of parameter names. In a one-compartment intravenous bolus model, dose enters the central compartment immediately, concentration is amount divided by volume, and elimination is first order. The corresponding nlmixr2 function needs to express those relationships, describe how parameters vary between individuals, and specify how observations differ from predictions.

This tutorial focuses on that implementation. Assume doses are in milligrams, time in hours, clearance in liters per hour, volume in liters, and observed concentrations in milligrams per liter. The model assumes complete instantaneous delivery of the dose into one well-mixed compartment, linear elimination, and constant clearance and volume during the observation period. It has no absorption process.

Definition: For individual $i$, let $A_i(t)$ be the amount of drug in the central compartment, $CL_i$ clearance, and $V_i$ volume. Between doses, the one-compartment model is $dA_i(t)/dt=-(CL_i/V_i)A_i(t)$, and the predicted concentration is $C_i(t)=A_i(t)/V_i$. An intravenous bolus of dose $D_i$ at time zero sets $A_i(0)=D_i$.

The solution after a single bolus is $C_i(t)=(D_i/V_i)\exp[-(CL_i/V_i)t]$. Clearance and volume therefore determine both the initial concentration, through $D_i/V_i$, and the rate of decline, through $CL_i/V_i$. The software implementation should encode amount in the compartment and calculate concentration from that amount; dose is not itself a concentration.

Worked Example: Checking a Structural Prediction

Suppose a hypothetical individual receives a 100 mg bolus, with clearance 5 L/h and volume 50 L. The elimination rate constant is $k=CL/V=5/50=0.1000$ per hour, and the initial concentration is $D/V=100/50=2.000$ mg/L.

At 4 hours, the model predicts:

C(4)=\frac{100\ \mathrm{mg}}{50\ \mathrm{L}} \exp\left[-\left(\frac{5\ \mathrm{L/h}}{50\ \mathrm{L}}\right)(4\ \mathrm{h})\right] =2.000\exp(-0.4000)=1.341\ \mathrm{mg/L}.

The exponent is dimensionless because the rate in per hour is multiplied by time in hours. The predicted concentration is lower than the initial concentration, as expected for positive clearance.

Typical Parameters and Interindividual Effects

A common population model represents individual clearance and volume as positive typical values multiplied by exponentiated random effects:

CL_i=\exp(\theta_{CL}+\eta_{CL,i}), \qquad V_i=\exp(\theta_V+\eta_{V,i}).

Here $\theta_{CL}$ and $\theta_V$ are the logarithms of the typical clearance and volume, respectively. The random effects $\eta_{CL,i}$ and $\eta_{V,i}$ describe individual deviations on the log scale. A standard assumption is that the random-effects vector has mean zero and a multivariate normal distribution. If the two effects are specified separately in the model function below, they have a diagonal covariance structure: each has its own variance and their covariance is set to zero. This is an assumption of the specification, not a finding established by the data.

In nlmixr2, assignments in the ini block give starting values for estimation. A declaration such as eta.cl ~ 0.1 supplies an initial value for the random-effect variance; it does not mean that the standard deviation is 0.1. The exponential link makes clearance and volume positive for all finite random-effect values.

Worked Example: Applying Individual Random Effects

For a hypothetical individual, let typical clearance be 5 L/h and typical volume be 50 L. Suppose $\eta_{CL}=\log(1.2)=0.1823$ and $\eta_V=\log(0.9)=-0.1054$. Then:

CL_i=5\exp(0.1823)\approx 6.000\ \mathrm{L/h}, \qquad V_i=50\exp(-0.1054)\approx 45.00\ \mathrm{L}.

For a 100 mg bolus, the predicted concentration at 2 hours is:

C_i(2)=\frac{100}{45.00}\exp\left[-\left(\frac{6.000}{45.00}\right)(2)\right] =2.222\exp(-0.2667)\approx 1.702\ \mathrm{mg/L}.

The concentration uses the individual’s clearance and volume, not the typical values. The initial concentration is $100/45=2.222$ mg/L, and the larger individual clearance produces a faster decline than would occur at the typical clearance, holding volume fixed.

Writing the Model Function

The following function states the model using a central amount, centr, and a predicted concentration, cp. The starting values are illustrative and are not claimed to be estimates from a fitted data set.

Model function:
one_cmt <- function() {
  ini({
    tcl <- log(5)
    tv <- log(50)
    eta.cl ~ 0.1
    eta.v ~ 0.1
    add.err <- 0.2
  })
  model({
    cl <- exp(tcl + eta.cl)
    v <- exp(tv + eta.v)
    d/dt(centr) = -cl / v * centr
    cp = centr / v
    cp ~ add(add.err)
  })
}

The initial typical values are $CL=\exp(\log 5)=5$ L/h and $V=\exp(\log 50)=50$ L. The differential equation describes amount declining by first-order elimination. The line defining cp converts that amount to concentration. The observation statement requests an additive residual error: conditional on the model prediction and individual effects, the observed concentration is modeled as the prediction plus an error with mean zero and standard deviation add.err. Its units are mg/L. In the initialization block, 0.2 mg/L is a starting value for that standard deviation, not a fixed known error.

The function does not create a bolus merely by defining the differential equation. The dose must be represented in the input data as a dose event assigned to the central compartment. In the example below, CMT=1 identifies that compartment, while EVID=1 identifies a dose record and EVID=0 identifies an observation record.

Worked Example: Preparing Bolus Records and Calling the Model

Suppose three hypothetical individuals each receive 100 mg at time zero, with concentration observations at 0.5, 1, 2, 4, and 8 hours. The following code shows the record structure. The concentrations are illustrative values; they are not the result of fitting the model.

Observation records:
obs <- data.frame(
  ID = rep(1:3, each = 5),
  TIME = rep(c(0.5, 1, 2, 4, 8), 3),
  DV = c(1.90, 1.82, 1.65, 1.42, 1.05,
       2.02, 1.88, 1.76, 1.49, 1.12,
       1.84, 1.75, 1.62, 1.37, 1.00),
  AMT = NA_real_, EVID = 0, CMT = 1
)

Dose records and combined data:
dose <- data.frame(ID = 1:3, TIME = 0, DV = NA_real_,
  AMT = 100, EVID = 1, CMT = 1)
dat <- rbind(dose, obs)
dat <- dat[order(dat$ID, dat$TIME, -dat$EVID), ]

Fit call:
fit <- nlmixr2(one_cmt, dat, est = "focei")

Each dose row has a positive amount, the dose event code, and no observed concentration. Each observation row has a measured concentration in DV, no dose amount, and the observation event code. Sorting puts each individual’s dose record before observations at later times; when records share a time, the chosen ordering places the dose before the observation. The fit call requests estimation with FOCEI. The ini values are starting values for that estimation, so they must not be reported as fitted results.

Before fitting, check that IDs are consistent, times are in hours, dose amounts are in milligrams, concentrations are in mg/L, and the dose is assigned to the compartment represented by centr. Also verify that dose records are not accidentally treated as concentration observations.

Residual Error and Basic Diagnostics

With additive error, the observation model is $DV_{ij}=cp_{ij}+\epsilon_{ij}$, where $j$ indexes observations and the residuals are assumed to have mean zero, constant variance, and to be independent conditional on the model and individual effects. This assumption means the same additive-error standard deviation applies at low and high predicted concentrations. It may be unsuitable when variability grows with concentration; a proportional or combined error form would then be a candidate to evaluate, not an automatic improvement.

For a simple residual check, compare observed values with predictions and examine residuals against time and prediction. A systematic decline or rise in residuals over time may indicate a structural mismatch, an incorrect dose time, or a data problem. A residual spread that grows with predicted concentration challenges constant additive variance. A few extreme residuals warrant checking records and units before changing the model. These checks diagnose possible problems; they do not by themselves identify a unique replacement model.

Worked Example: Interpreting an Additive Residual

Suppose the model predicts 1.702 mg/L for one observation and the hypothetical measured value is 1.782 mg/L. The raw residual, defined here as observed minus predicted, is:

1.782-1.702=0.080\ \mathrm{mg/L}.

If the illustrative additive-error standard deviation is 0.200 mg/L, the residual divided by that standard deviation is $0.080/0.200=0.400$. This is a modest positive discrepancy relative to the assumed residual scale. It is not evidence by itself that the model is correct; interpretation requires examining residual patterns across the data and the fitted model.

Derivation, Failure Modes, and Limits

The model’s differential equation also verifies that the concentration prediction used above is consistent with the implemented state. With constant $CL_i$ and $V_i$, divide the amount equation by $V_i$:

\frac{dC_i(t)}{dt} =\frac{1}{V_i}\frac{dA_i(t)}{dt} =-\frac{CL_i}{V_i}\frac{A_i(t)}{V_i} \quad\text{with}\quad C_i(t)=\frac{A_i(t)}{V_i}.

Thus $dC_i(t)/dt=-(CL_i/V_i)C_i(t)$, whose solution with $C_i(0)=D_i/V_i$ is $C_i(t)=(D_i/V_i)\exp[-(CL_i/V_i)t]$. This establishes why the amount equation, the concentration calculation, and the worked predictions describe the same one-compartment model.

A frequent implementation mistake is to put a dose into a model state that has not been matched to the dose record, then conclude that the model predicts near-zero concentrations because clearance or volume is wrong. The correction is to check event coding and compartment assignment first: the bolus must enter centr, the amount state used in the differential equation. Another mistake is to treat $CL/V$ as clearance. It is an elimination rate with units of inverse time; clearance has units of volume per time.

This implementation also has limits. It assumes instantaneous intravenous input, one well-mixed compartment, linear elimination, and the specified random-effect and residual-error distributions. It does not represent absorption, distribution into additional compartments, time-varying parameters, or a different observation process. If these assumptions are implausible, fitting the function successfully does not make its assumptions true. The additive model is one explicit residual choice; its constant absolute error scale should be judged against data patterns and scientific context.

Key takeaway: In this nlmixr2 implementation, the dose record supplies the amount to the central state, clearance and volume determine its decline, $cp=centr/v$ defines the prediction, exponential random effects preserve positive individual parameters, and the additive statement specifies residual variation in concentration units.

Check Your Understanding

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

For a 60 mg bolus, clearance 3 L/h, and volume 30 L, calculate the predicted concentration at 2 hours.
If $tcl=\log(4)$ and $\eta_{CL}=\log(1.25)$, what individual clearance does the model imply?
In the example data structure, what are the key fields for distinguishing a dose record from an observation record?
What does an initial value of 0.1 in eta.cl ~ 0.1 represent, and is it a standard deviation?
Name one diagnostic pattern that would make an additive residual-error assumption questionable.

Solutions. (1) $k=3/30=0.1000$ per hour, so $C(2)=(60/30)\exp[-0.1000(2)]=2\exp(-0.2)=1.637$ mg/L. (2) $CL_i=\exp[\log(4)+\log(1.25)]=4(1.25)=5.000$ L/h. (3) A dose row has EVID=1, an amount in AMT, and no DV; an observation row has EVID=0 and a measured value in DV. Both are assigned to the appropriate compartment. (4) It is an initial value for the random-effect variance, not a standard deviation; the corresponding initial standard deviation would be $\sqrt{0.1}$. (5) Residual spread that increases with predicted concentration challenges the constant absolute variance assumed by additive error.