Why a Unit Error Can Look Like a Physiological Failure
A PBPK simulation can produce a curve that looks dramatically wrong even when its equations have the intended physiological structure. A concentration that falls almost immediately, an organ that accumulates drug far too slowly, or an amount that becomes negative may reflect a mismatch between the units used for time, volume, or rate rather than an unusual biological process. Such errors are especially easy to introduce when parameter values and model code come from sources that use different unit conventions.
The earlier material on dimensional analysis and amount-to-concentration conversion supplies the basic checks: every term in an equation must have compatible units, and an amount divided by a volume is a concentration. Here the emphasis is diagnostic. Start from the implausible output, identify the equation that controls it, and compare the units actually used by the code with the units required by that equation.
A useful first step is to write down the unit of each quantity as implemented, not merely as intended. In the examples below, \(A\) denotes amount in mg, \(C\) denotes concentration in mg/L, \(V\) denotes volume in L, \(Q\) denotes flow in L/h, \(CL\) denotes clearance in L/h, and time is expressed in hours unless stated otherwise. The same diagnostic logic applies when another coherent set of units is used, provided conversions are made consistently.
Trace the Scale of a Suspicious Change
A discrepancy by a familiar conversion factor is a strong clue. A change that is about 60-fold too fast or slow suggests hours-versus-minutes; a change about 1,000-fold off often points to liters-versus-milliliters. These clues are not proof by themselves: parameters may also be wrong, and models can be physiologically misspecified. They do, however, help prioritize which unit pathway to inspect.
For a well-mixed compartment with amount \(A\), volume \(V\), and clearance \(CL\), elimination is commonly written \(dA/dt=-CL\,C\), where \(C=A/V\). With constant \(CL\) and \(V\), this gives \(dA/dt=-(CL/V)A\). Thus the elimination rate constant \(CL/V\) must have units of inverse time. A mismatch between the time unit for \(CL\) and the time unit used by the solver directly changes the simulated rate.
Worked Example: A Concentration Falls Too Quickly
Model and assumptions. Consider a hypothetical one-compartment model with instantaneous mixing, constant volume, and first-order elimination. Let \(V=12\) L and \(CL=1.2\) L/h. There is no ongoing input. The intended time unit is hours.
The elimination rate constant is
Suppose the simulation clock advances in minutes, but the clearance value is entered as the unchanged number 1.2 and is therefore treated as 1.2 L/min. Then the implemented rate constant is
After 60 minutes, the intended remaining fraction is \(e^{-0.100\ \mathrm{h}^{-1}\times1\ \mathrm{h}}\), whereas the incorrectly implemented remaining fraction is \(e^{-0.100\ \mathrm{min}^{-1}\times60\ \mathrm{min}}\). Numerically,
For a hypothetical initial amount of 120 mg, these correspond to 108.6 mg and 0.2975 mg remaining, respectively. The model has eliminated drug 60 times faster than intended because the numerical clearance was interpreted per minute rather than per hour. If the solver uses minutes, convert the clearance: \(1.2\ \mathrm{L/h}=0.0200\ \mathrm{L/min}\). The resulting rate constant is \(0.001667\ \mathrm{min}^{-1}\), and \(e^{-0.001667\times60}=0.9048\), as required.
This calculation also demonstrates a quick diagnostic: calculate the implied rate constant from the implemented values and check its unit. If a parameter is documented per hour but the implementation treats it as per minute, the code may still run without warning, because its arithmetic is numerically possible. The output, however, represents a different model.
Find Volume-Scale Errors in Flow-Driven Exchange
For a well-mixed compartment receiving flow \(Q\), a common exchange term for concentration is \(Q(C_{\mathrm{in}}-C)\). Since \(Q\) has units L/h and concentration has units mg/L, the product has units mg/h, appropriate for an amount rate. The corresponding concentration equation divides the amount rate by compartment volume: \(dC/dt=(Q/V)(C_{\mathrm{in}}-C)\). The ratio \(Q/V\) must be in inverse time, with numerator and denominator expressed using compatible volume units.
Worked Example: A Compartment Exchanges Drug Too Slowly
Model and assumptions. Consider a hypothetical, well-mixed compartment of volume \(V=2.00\) L supplied by constant flow \(Q=5.00\) L/h. Assume the inlet concentration \(C_{\mathrm{in}}\) is constant, there is no elimination within the compartment, and the initial concentration is zero. The exchange equation is \(dC/dt=(Q/V)(C_{\mathrm{in}}-C)\).
With consistent units, \(Q/V=5.00/2.00=2.50\ \mathrm{h}^{-1}\). Its solution from \(C(0)=0\) is \(C(t)=C_{\mathrm{in}}(1-e^{-(Q/V)t})\). After one hour, the fraction of inlet concentration reached is \(1-e^{-2.50}=0.9179\).
Now suppose the volume is entered as 2000, because it was recorded in mL, while the flow remains numerically 5.00 in L/h. The code forms \(5.00/2000=0.00250\), and may treat this as \(0.00250\ \mathrm{h}^{-1}\). It then predicts a one-hour fraction \(1-e^{-0.00250}=0.002497\), or about 0.250% of the inlet concentration, rather than 91.8%. The numerical ratio is 1,000 times smaller because the volume number was entered in milliliters while the flow number was in liters per hour.
The dimensional check reveals the defect: \(5.00\ \mathrm{L/h}\) divided by \(2000\ \mathrm{mL}\) is not a rate in inverse hours. Convert the volume to 2.00 L before using the ratio, or convert the flow to 5000 mL/h and use the volume 2000 mL. Either consistent choice gives \(2.50\ \mathrm{h}^{-1}\). This check also guards against the opposite error, in which a flow is converted to mL/h but a volume remains in L, creating a rate that is 1,000 times too large.
When the volume error affects an amount-to-concentration conversion rather than a flow-to-volume ratio, the symptom can be a concentration that is 1,000 times too small or too large. For example, 100 mg in 2 L is 50 mg/L; in 2000 mL it is 0.050 mg/mL. Those are the same concentration expressed in different units. Treating 0.050 mg/mL as though it were 0.050 mg/L introduces a 1,000-fold reporting error.
Check Whether the State Is an Amount or a Concentration
Rates of input must match the state they change. If the state is amount \(A\), an infusion rate \(R\) in mg/h can enter \(dA/dt\) directly. If the state is concentration \(C\) in mg/L, the same input must be divided by the relevant volume: \(dC/dt=R/V\), with units \((\mathrm{mg/h})/\mathrm{L}=\mathrm{mg/(L\,h)}\). Adding a rate in mg/h directly to a concentration derivative in mg/(L h) is dimensionally invalid.
Worked Example: An Amount Input Is Added to a Concentration State
Model and assumptions. Consider a hypothetical 3.00 L, well-mixed compartment whose state variable is concentration \(C\) in mg/L. A constant infusion supplies \(R=3.00\) mg/h for two hours. Assume no elimination, no exchange with other compartments, and an initial concentration of zero.
The input must be divided by volume. Therefore,
After two hours, \(C=0+(1.00\ \mathrm{mg/(L\,h)})(2.00\ \mathrm{h})=2.00\) mg/L. Equivalently, the infusion delivers 6.00 mg, and \(6.00\ \mathrm{mg}/3.00\ \mathrm{L}=2.00\) mg/L.
If the code instead adds the numerical value 3.00 directly to \(dC/dt\), it implicitly treats the input as 3.00 mg/(L h), although the supplied parameter is 3.00 mg/h. The simulated increase is then \((3.00\ \mathrm{mg/(L\,h)})(2.00\ \mathrm{h})=6.00\) mg/L, three times the correct increase for this 3.00 L compartment. This particular factor happens to equal the volume in liters; the error will vary with volume and cannot be repaired by a universal multiplier.
Negative Amounts: Separate Unit Errors From Numerical Overshoot
A physically nonnegative amount becoming negative is a warning, not a diagnosis. It may result from a wrongly scaled elimination rate, an inconsistent time step, a sign or flux-bookkeeping error, or numerical integration error. For the exact first-order elimination equation \(dA/dt=-kA\), with \(A(0)>0\) and \(k\geq0\), the solution is \(A(t)=A(0)e^{-kt}\), which stays positive at every finite time. A negative computed amount therefore cannot be the exact solution of this model under those assumptions.
A simple derivation shows how a numerical method can nevertheless produce a negative value. Explicit Euler approximates the solution by \(A_{j+1}=A_j-kA_j\Delta t=A_j(1-k\Delta t)\). If \(A_j>0\), the next value is nonnegative only when \(0\leq k\Delta t\leq1\). A rate expressed in inverse minutes combined with a step interpreted as 60 minutes can make \(k\Delta t\) large, causing the update to cross zero. This condition describes explicit Euler for this equation; it is not a universal guarantee for every solver or every PBPK system.
Worked Example: Distinguishing a Negative Euler Update From the Exact Model
Model and assumptions. Take a hypothetical amount \(A_0=100\) mg undergoing first-order elimination with \(k=0.100\ \mathrm{min}^{-1}\). Assume there is no input. Consider one explicit Euler step of \(\Delta t=60\) minutes.
The dimensionless step product is \(k\Delta t=0.100\times60=6.00\). Euler gives \(A_1=100(1-6.00)=-500\) mg. By contrast, the exact solution after 60 minutes is \(A(60)=100e^{-6}=0.2479\) mg, which is positive. The large negative Euler result is an integration overshoot for this step size, not a physically meaningful prediction.
If the intended parameter was instead \(0.100\ \mathrm{h}^{-1}\), a 60-minute interval is one hour, so \(k\Delta t=0.100\) when time is expressed in hours. The Euler update is \(100(1-0.100)=90.0\) mg, while the exact one-hour amount is \(100e^{-0.100}=90.48\) mg. This is a modest approximation error, not a negative amount. The calculation shows why both time-unit consistency and solver-step adequacy should be checked.
A Practical Diagnostic Sequence
Use the following sequence to localize the problem before changing physiological parameters. It combines equation-level checks with checks of the actual simulated trajectory.
Identify which amount or concentration first departs from expectation, and inspect the equation that updates it. This narrows the search to relevant inputs, flows, clearances, volumes, and exchange terms.
Record the units used in code for each state, parameter, and time variable. Do not rely solely on a parameter file's label; check how the numerical value is passed into the equation and solver.
An amount derivative must have units of amount per time; a concentration derivative must have units of concentration per time. Check products such as flow times concentration and ratios such as clearance divided by volume.
Estimate a rate constant, a concentration from amount divided by volume, or total infused amount over the interval. Look for factors near 60 or 1,000 and compare the model with a simple case whose expected behavior is clear.
Repeat the run with tighter tolerances or a smaller maximum step, and inspect solver warnings and mass-balance residuals. If results change materially with these settings, numerical integration needs attention as well as unit consistency.
What the Diagnostics Can and Cannot Establish
Dimensional analysis is decisive when an equation does not have the units required for its state derivative. It is less decisive when code uses dimensionless numbers with an implicit convention: the code may be internally consistent but represent minutes where hours were intended. In that case, parameter documentation, time-grid definitions, input schedules, and calculated rate scales must be compared. A simulation that finishes without an error message is not evidence that its units are correct.
Unit checks complement, but do not replace, other debugging methods. Mass-balance residuals can reveal unexplained gains or losses of drug, but a model may conserve mass while using a clearance that is 60 times too large. Physiological plausibility checks can flag a curve as suspicious, but cannot identify which conversion caused it. Solver diagnostics can reveal sensitivity to step size or tolerance, but tightening tolerances will not correct a clearance entered in the wrong time unit. Use these approaches together rather than treating any one as conclusive.
A frequent mistake is to see rapid elimination, reduce clearance until the curve looks reasonable, and stop investigating. That adjustment can hide an hours-versus-minutes error and leave other terms—such as infusion schedules and flows—on inconsistent time scales. The correction is to verify the units and conversions first, then assess whether the remaining simulation is physiologically and numerically credible. Likewise, do not interpret a negative amount as proof of a unit error alone: inspect signs, conservation, step size, and state definitions.
Check Your Understanding
For each question, identify the main unit or numerical issue before reading the solution.
- A model uses \(CL=0.60\) L/h and \(V=6.0\) L, but advances time in minutes while leaving the clearance value unchanged. What rate constant does the code imply, and how does it compare with the intended rate?
- A concentration state is in mg/L, and an input of 12 mg/h is applied to a 4 L compartment. What input term belongs in its concentration derivative?
- A compartment has \(Q=2\) L/h and \(V=500\) mL. What is the correct value of \(Q/V\) in inverse hours?
- For explicit Euler applied to \(dA/dt=-kA\), what condition on \(k\Delta t\) keeps the next amount nonnegative when the current amount is positive? Does failure of this condition prove that the model's exact amount is negative?
Solutions. (1) The intended rate constant is \(0.60/6.0=0.100\ \mathrm{h}^{-1}\). Unconverted, the code interprets it as \(0.100\ \mathrm{min}^{-1}\), 60 times the intended rate per hour; the clearance should be \(0.0100\) L/min for a minute-based clock. (2) Divide the amount rate by volume: \(12\ \mathrm{mg/h}/4\ \mathrm{L}=3.00\ \mathrm{mg/(L\,h)}\). (3) Convert 500 mL to 0.500 L, giving \(2/0.500=4.00\ \mathrm{h}^{-1}\). (4) Euler gives \(A_{\mathrm{next}}=A(1-k\Delta t)\), so \(0\leq k\Delta t\leq1\) is required for a nonnegative next value. Violating this condition can make the Euler update negative, but the exact solution \(A(t)=A(0)e^{-kt}\) remains positive for positive \(A(0)\) and nonnegative \(k\).