Checking Simulated Concentration-Time Profiles for Model Consistency

Why Check the Profile, Not Only the Exposure Summary?

A simulated exposure summary can look plausible even when the concentration-time curve that produced it does not follow the specified model. A concentration may be expressed in the correct units, and an area under the curve may be within an expected range, while a dose has been applied at the wrong time, an absorption rate has been mis-coded, or a model component has been omitted. Inspecting profiles is therefore a direct check of whether the simulation behaves as its structural model says it should.

The check is conditional on the model and regimen being examined. It is not a search for one universally “correct” curve: an oral dose with first-order absorption, an intravenous bolus in a one-compartment model, and an intravenous bolus in a two-compartment model imply different shapes. Before judging a simulated profile, specify the route, dose time, units, parameter values, and whether the profile is a latent model prediction or includes simulated residual error.

Definition: A model-consistency check compares a simulated concentration-time profile with the qualitative shape and quantitative relationships implied by the stated pharmacokinetic model, parameter values, route, and dose history. A mismatch is evidence to investigate, not by itself proof of a particular coding error.

What Shape Does the Model Imply?

For an intravenous bolus in a one-compartment model with linear elimination, concentration begins at dose divided by apparent volume and declines exponentially. If \(D\) is dose, \(V\) is volume, and \(k_e=CL/V\) is the elimination rate constant, then after a single bolus at time zero:

C(t)=\frac{D}{V}e^{-k_e t}, \qquad t\geq 0.

This curve has an immediate post-dose concentration \(D/V\), is positive, and decreases continuously when \(k_e>0\). Taking logarithms gives \(\ln C(t)=\ln(D/V)-k_e t\): on a plot of log concentration against time, the model predicts a straight line with slope \(-k_e\). This straight-line expectation applies to this model and this post-bolus phase; it is not a rule for every route or compartment structure.

For a one-compartment model with first-order absorption after an oral dose, assuming a fixed lag-free dose at time zero, bioavailability \(F\), first-order absorption rate \(k_a\), and first-order elimination rate \(k_e\), the concentration is:

C(t)=\frac{F D k_a}{V(k_a-k_e)} \left(e^{-k_e t}-e^{-k_a t}\right), \qquad t\geq 0,\quad k_a\neq k_e.

The curve starts at zero, rises while drug is entering the systemic circulation faster than it is being eliminated, reaches a peak, and then declines. The peak time follows by differentiating this expression and setting the derivative to zero. The derivative is proportional to \(-k_e e^{-k_e t}+k_a e^{-k_a t}\). Thus, at the peak, \(k_a e^{-k_a t}=k_e e^{-k_e t}\), which gives:

t_{\max}=\frac{\ln(k_a/k_e)}{k_a-k_e}.

For positive, unequal \(k_a\) and \(k_e\), this time is positive. If \(k_a>k_e\), absorption is faster than elimination and the terminal decline is governed by \(k_e\). If \(k_a<k_e\), the slower absorption can govern the apparent terminal decline, a phenomenon called flip-flop kinetics. Therefore, a terminal slope that differs from the elimination rate is not automatically an implementation error; the structural model and rate ordering matter.

A model with distribution between a central and a peripheral compartment can produce a rapid early decline followed by a slower terminal decline after an intravenous bolus. A common representation is a sum of exponentials, \(C(t)=A e^{-\alpha t}+B e^{-\beta t}\), where \(A,B>0\) and \(\alpha>\beta>0\). Early concentration changes can be dominated by the \(\alpha\) term, while sufficiently late concentrations are dominated by the \(\beta\) term. A single straight line on the entire log-concentration plot is not expected in this case.

A Practical Profile-Checking Sequence

Confirm the simulation setup.
Check dose amount and units, route, dose and sampling times, parameter units, and the model used to generate the profile. Confirm whether plotted values are latent concentrations or include residual error.

Check basic boundary and sign behavior.
For a lag-free oral first-order absorption profile, concentration should start at zero. For an intravenous bolus, the immediate post-dose concentration should agree with dose divided by volume in the one-compartment model. Concentrations should be nonnegative.

Compare timing and phases.
For oral dosing, check whether the rise, peak, and decline agree with the absorption and elimination rates. For multi-compartment dosing, look for the expected early and late phases rather than forcing one slope across the full curve.

Quantify discrepancies.
Calculate model-predicted concentrations at selected times, or compare log concentration slopes over intervals where the model predicts a single exponential phase. Make comparisons at the same times and under the same dosing conditions.

Trace the mismatch before changing the model.
First investigate time-unit conversion, dose timing, parameter mapping, output columns, and numerical implementation. Change model structure only when the mismatch is systematic and the intended structural assumptions themselves are inadequate.

For simulated data with residual error, an individual observed profile need not lie exactly on its latent model curve. Check the latent prediction separately when possible, and assess noisy simulated observations against the error model rather than treating every fluctuation as structural failure. The time grid also matters: a coarse grid can miss a sharp peak without the underlying model being wrong.

Worked Example: Checking an Oral Absorption Profile

Consider a hypothetical, single oral dose whose concentration model is \(C(t)=12(e^{-0.2t}-e^{-1.2t})\) mg/L, with time \(t\) in hours. This form corresponds to a lag-free, first-order absorption and elimination model with \(k_e=0.2\) h\(^{-1}\), \(k_a=1.2\) h\(^{-1}\), and a scale factor of 12 mg/L. The rates are positive and unequal, and \(k_a>k_e\), so the curve should rise and then have a terminal decline governed by elimination.

At \(t=1\) h, the predicted concentration is:

C(1)=12(e^{-0.2}-e^{-1.2}) =12(0.81873075-0.30119421) =6.2104385\ \text{mg/L} \approx 6.210\ \text{mg/L}.

The expected peak time is:

t_{\max}=\frac{\ln(1.2/0.2)}{1.2-0.2} =\frac{\ln 6}{1.0} =1.7918\ \text{h}.

A simulated latent profile that is already declining from its first post-dose value, or whose peak is far from 1.792 h, deserves investigation. Possible causes include swapped rate parameters, a time recorded in minutes but interpreted as hours, or an incorrectly applied dose time. The calculation does not identify which cause occurred; it supplies a model-based reference for the check.

Checking Exponential Decline and Distribution Phases

For a one-compartment intravenous bolus profile, an interval-specific elimination rate can be calculated from two positive concentrations \(C(t_1)\) and \(C(t_2)\), with \(t_2>t_1\):

\widehat{k}_{e,\mathrm{interval}} =\frac{\ln C(t_1)-\ln C(t_2)}{t_2-t_1}.

This is simply the negative slope of the line joining the two log concentrations. Under the stated one-compartment model, without ongoing input, it should equal \(k_e\) apart from rounding or any numerical approximation. With noisy observations it is an interval estimate affected by residual error; with a multi-compartment profile it can vary by interval because the relative contributions of the exponential phases change.

Worked Example: Detecting a Late-Interval Decline Mismatch

A hypothetical one-compartment intravenous bolus has \(D=100\) mg, \(V=20\) L, and \(k_e=0.15\) h\(^{-1}\). Assume a single dose at time zero, no ongoing input, and concentrations in the simulation are latent model predictions. The initial concentration should be \(D/V=100/20=5.00\) mg/L. The model predicts:

C(4)=5e^{-0.15(4)}=5e^{-0.6} =5(0.5488116)=2.744\ \text{mg/L},

and:

C(8)=5e^{-0.15(8)}=5e^{-1.2} =5(0.3011942)=1.506\ \text{mg/L}.

Suppose the simulated values are 5.00 mg/L at time zero, 2.74 mg/L at 4 h, and 2.40 mg/L at 8 h. The initial value and the first interval agree closely with the model, but the 4-to-8 h decline does not. Its implied interval rate is:

\widehat{k}_{e,\mathrm{interval}} =\frac{\ln(2.74)-\ln(2.40)}{8-4} =\frac{\ln(2.74/2.40)}{4} \approx 0.03312\ \text{h}^{-1}.

This is much smaller than the specified \(0.15\) h\(^{-1}\). For this deterministic profile, that discrepancy is inconsistent with the stated one-compartment, post-bolus model. Check whether the elimination parameter was passed correctly, whether time units were converted consistently, and whether the displayed values came from the intended simulation. Do not infer from the discrepancy alone that the biological model needs a second compartment; implementation and output errors are plausible alternatives.

Worked Example: Recognizing a Two-Phase Decline

Consider a hypothetical central concentration after an intravenous bolus in a two-compartment model, represented here as \(C(t)=4e^{-1.0t}+1e^{-0.1t}\) mg/L, with time in hours. Both exponential coefficients are positive and the faster rate, \(1.0\) h\(^{-1}\), exceeds the slower rate, \(0.1\) h\(^{-1}\). The initial concentration is \(4+1=5\) mg/L.

At 4 h:

C(4)=4e^{-4}+e^{-0.4} =4(0.0183156)+0.6703200 =0.7435824\ \text{mg/L}.

At 8 h:

C(8)=4e^{-8}+e^{-0.8} =4(0.00033546)+0.4493290 =0.4506708\ \text{mg/L}.

At 12 h:

C(12)=4e^{-12}+e^{-1.2} =4(0.000006144)+0.3011942 =0.3012188\ \text{mg/L}.

The 4-to-8 h interval rate is \(\ln(0.7435824/0.4506708)/4\approx0.125\) h\(^{-1}\), whereas the 8-to-12 h interval rate is \(\ln(0.4506708/0.3012188)/4\approx0.1007\) h\(^{-1}\). The later estimate is close to the slower rate of \(0.1\) h\(^{-1}\), as expected once the fast component is small; the interval estimate is approximately \(0.1007\) h\(^{-1}\). A curved log profile that becomes approximately linear late is therefore compatible with this two-phase model; demanding a single slope over all times would misdiagnose the model.

Diagnostics, Alternatives, and Limits

A useful diagnostic set pairs the concentration-time plot on its original scale with a log concentration-time plot when concentrations are positive. The original scale makes peak timing and absolute concentrations visible; the log scale makes exponential phases easier to inspect. Tabulated predictions at selected times can help distinguish a visual impression from a numerical mismatch. For oral dosing, include points before and around the predicted peak. For an expected terminal phase, include sufficiently late points, while recognizing that concentrations near numerical or reporting limits may be unreliable.

Common patterns point to different checks, but none uniquely identifies the cause:

Incorrect initial value after an intravenous bolus: verify dose units, volume units, and whether the plotted time-zero value means immediately before or immediately after the dose.
Oral profile begins above zero or peaks at an implausible time: inspect lag-time handling, absorption input, and whether \(k_a\) and \(k_e\) were assigned to the intended parameters.
Decline is too shallow or too steep: verify time units and clearance/volume parameterization before interpreting the slope as evidence about model structure.
Unexpected oscillations or negative concentrations: check numerical implementation, interpolation, output transformation, and the model equations. A smooth positive latent concentration is expected for the positive-parameter examples above.
One subject differs while others behave as expected: inspect that subject’s sampled parameters and dose history, and distinguish a valid individual parameter combination from a programming error.

Direct comparison with the structural equations is the strongest check when the model and parameters are known. Shape-based inspection is quicker and can reveal timing errors, but it is less specific. Checking derived features such as \(t_{\max}\), initial concentration, or interval slopes is more quantitative, although it can miss errors between sampled times and can be sensitive to noise. If the simulation includes residual error, compare its observations with the distribution implied by the error model rather than expecting exact equality to the latent curve. The main limitation is identifiability: several different implementation faults can create similar profile discrepancies, so follow-up checks of code, units, and parameter inputs are needed.

Common mistake: “The log concentration plot is not one straight line, so the simulation must be wrong.” This conclusion is valid only when the specified model predicts a single exponential decline over the inspected interval. A two-compartment profile can have an early distribution phase and a later elimination phase; oral absorption can also shape the apparent terminal slope. First state the model’s predicted phases, then compare the relevant part of the profile with that prediction.

Key takeaway: A profile check is a conditional test of internal consistency: given the stated route, dose history, model structure, parameter values, and time units, do simulated concentrations have the expected boundaries, timing, phases, and rates? Diagnose setup and implementation before attributing every mismatch to model structure.

Check Your Understanding

Answer each question, then compare with the solutions below.

For the oral model in the worked example, what is the expected peak time, and why should the profile initially rise?
A one-compartment intravenous bolus model has \(D=60\) mg and \(V=12\) L. What is the immediate post-dose concentration?
A profile is generated with two positive exponential components having rates \(\alpha=0.8\) h\(^{-1}\) and \(\beta=0.08\) h\(^{-1}\). Which component is expected to dominate sufficiently late?
Why is a log concentration plot that curves early and becomes nearly straight late not necessarily evidence of a simulation error?
For the one-compartment model with \(k_e=0.2\) h\(^{-1}\), what concentration ratio is predicted over a 5 h interval with no ongoing input?

Solutions. (1) \(t_{\max}=\ln(1.2/0.2)/(1.2-0.2)=\ln(6)=1.792\) h. Initially the absorption input exceeds elimination, so concentration rises from zero. (2) \(D/V=60/12=5.00\) mg/L. (3) The slower \(\beta\) component dominates late because the faster exponential decays more quickly. (4) A multi-compartment model can produce a rapid early distribution phase followed by a slower terminal phase, so one straight line is not expected across all times. (5) \(C(t+5)/C(t)=e^{-k_e(5)}=e^{-1}=0.3679\), under the stated single-exponential, no-input conditions.