Checking Structural Equations Against Model Code

Why Derive the Equations Independently?

Reading code and its written specification side by side is useful, but it does not by itself establish that the implemented equations have the intended mathematical behavior. A reviewer who independently derives an expected amount or concentration equation has a reference that does not depend on the code under review. The comparison can reveal a missing input term, a sign error, a mistaken compartment, or a formula that has the right general shape but the wrong time dependence.

The goal is not to demand that code use the same notation or algebraic expression as a derivation. A differential-equation solver and a closed-form expression, for example, can represent the same model. The question is whether the code implies the same structural relationships under the same assumptions and conditions.

Definition: A structural equation check independently states the model’s differential equations, initial conditions, and relevant inputs, derives an expected amount or concentration, and compares that result with the corresponding quantity implied by the implemented code. The check concerns the model’s mathematical structure; it does not alone verify the data, estimation procedure, or every aspect of a software run.

In the examples below, \(A\) denotes an amount, \(C\) a concentration, \(V\) a volume, \(D\) a dose amount, and \(t\) time. A first-order rate constant is denoted by \(k\), \(k_a\), or \(k_e\), as appropriate. The examples use hypothetical values. Each equation is conditional on its stated model assumptions; it is not a claim that a particular model is appropriate for every drug or data set.

A Focused Comparison Process

Start by translating the intended model into a small, explicit mathematical statement. Identify the state variables, their initial values, the inputs, and the rates that change each state. Then derive a result that can be compared with the code, such as a closed-form equation or a derivative at a particular time. A single matching value is weak evidence: different equations can coincide at one time. Compare multiple times or a property such as the initial slope or long-time behavior as well.

Write the independent model.
State the differential equation for each amount, the initial conditions, and any input that changes over time.

Derive a reference quantity.
Solve the equations when a closed form is available, or derive a checkable property such as an initial derivative or mass-balance relationship.

Trace the code’s mathematical meaning.
Identify which code expressions define the state change and which calculate the reported amount or concentration. Do not infer equivalence from variable names alone.

Compare beyond one point.
Evaluate selected times and inspect initial conditions, slopes, and limiting behavior. Investigate any disagreement before treating the code output as the reference.

These checks assume that the mathematical interpretation of the code is clear, that the comparison uses the same parameter values and initial conditions, and that the model conditions (for example, whether an input is on or off) are aligned. If code uses a numerical solver, small discrepancies can arise from numerical tolerances and time stepping; a large or systematic discrepancy calls for investigation rather than being dismissed as numerical error.

One-Compartment IV Bolus: Deriving the Exponential Decline

Consider a one-compartment model after an IV bolus. Assume the dose is instantaneously placed in the compartment at time zero, elimination is first-order with a constant rate \(k\), and there are no additional inputs. The amount immediately after dosing is \(A(0)=D\). Concentration is calculated as \(C(t)=A(t)/V\), with constant \(V\).

Model: For \(t\geq 0\), \(\frac{dA}{dt}=-kA\), with \(A(0)=D\). Under these assumptions, the amount decreases in proportion to the amount present.

Separate variables and integrate from time zero to \(t\). Since \(A\) is positive for a positive dose, the logarithm is defined:

\frac{dA}{A}=-k\,dt,\qquad \int_D^{A(t)}\frac{dA}{A}=-k\int_0^t dt, \qquad \log\!\left(\frac{A(t)}{D}\right)=-kt.

Exponentiating gives \(A(t)=D e^{-kt}\), and therefore \(C(t)=D e^{-kt}/V\). This derivation makes two direct code checks possible: the derivative at time zero must be \(-kD\), and the amount must approach zero as \(t\) increases without bound when \(k>0\).

Worked Example: Checking an IV Bolus Equation

Suppose a hypothetical model has \(D=100\) mg, \(k=0.2500\) h\(^{-1}\), and \(V=20.0\) L. Its intended structural equation is \(dA/dt=-kA\), with \(A(0)=100\) mg. The assumptions are a single bolus, constant first-order elimination, and no further input.

The independent solution is \(A(t)=100e^{-0.2500t}\) mg. At \(t=4.00\) h:

A(4)=100e^{-0.2500(4)} =100e^{-1} =36.78794412\ldots\ \text{mg}.

Thus \(C(4)=A(4)/V=36.78794412\ldots/20.0=1.839397206\ldots\) mg/L, or 1.839 mg/L to four significant figures. At time zero, \(C(0)=100/20.0=5.000\) mg/L. The initial amount slope is \(-kA(0)=-(0.2500)(100)=-25.00\) mg/h.

An implementation whose state equation is \(dA/dt=-kA\) and whose reported concentration is \(A/V\) has the expected form. If the code instead uses \(dA/dt=+kA\), its initial slope is positive and its amount grows, contradicting the independently derived decline. If the derivative is correct but concentration is calculated from a different quantity, the state equation alone will not catch the reporting mismatch.

Constant-Rate Infusion: Check the Input and the Off-Infusion Phase

A constant-rate infusion adds a nonzero input to the amount equation while the infusion is running. Assume a one-compartment model, a constant infusion rate \(R_0\) from time zero through duration \(T\), constant first-order elimination \(k\), and \(A(0)=0\). The input is \(R_0\) during infusion and zero after it stops. This example treats the switch at \(T\) as an imposed model condition; it does not address how dose events are represented in a particular software system.

Model: During infusion, \(\frac{dA}{dt}=R_0-kA\). After infusion stops, \(\frac{dA}{dt}=-kA\). The infusion-phase solution is \(A(t)=\frac{R_0}{k}(1-e^{-kt})\) for \(0\leq t\leq T\); after \(T\), \(A(t)=A(T)e^{-k(t-T)}\).

For the infusion phase, multiply the differential equation by the integrating factor \(e^{kt}\). The left side becomes the derivative of \(e^{kt}A(t)\). Integrating and using \(A(0)=0\) gives the result:

\frac{d}{dt}\!\left(e^{kt}A(t)\right)=R_0e^{kt}, \qquad e^{kt}A(t)-A(0)=\frac{R_0}{k}(e^{kt}-1), \qquad A(t)=\frac{R_0}{k}(1-e^{-kt}).

Once the input is zero, the amount already present declines according to the no-input equation. The two expressions agree at \(t=T\): the post-infusion expression gives \(A(T)e^{-k(T-T)}=A(T)\). A code check that examines only the infusion phase could miss an implementation that continues the input after the intended stopping time.

Worked Example: Checking Both Phases of an Infusion

For a hypothetical infusion, let \(R_0=30.0\) mg/h, \(k=0.2000\) h\(^{-1}\), \(T=3.00\) h, and \(V=15.0\) L. Assume the amount starts at zero and the input stops at \(T\). The constant-rate input and first-order elimination assumptions justify the two equations above.

At \(t=2.00\) h, infusion is running:

A(2)=\frac{30.0}{0.2000}(1-e^{-0.4000}) =150(1-0.6703200460\ldots) =49.45199310\ldots\ \text{mg}.

The corresponding concentration is \(49.45199310\ldots/15.0=3.296799540\ldots\) mg/L, or 3.297 mg/L. At the stopping time, \(A(3)=150(1-e^{-0.6000})=67.67825459\ldots\) mg. At \(t=5.00\) h, two hours after stopping:

A(5)=A(3)e^{-0.2000(5-3)} =150(e^{-0.4000}-e^{-1.0000}) =45.36609072\ldots\ \text{mg}.

So \(C(5)=45.36609072\ldots/15.0=3.024406048\ldots\) mg/L, or 3.024 mg/L. The amount is continuous at \(T\), while its derivative changes because the input changes from \(R_0\) to zero. For example, just before stopping the derivative is \(30.0-(0.2000)A(3)\), approximately \(16.46\) mg/h; just after stopping it is \(-(0.2000)A(3)\), approximately \(-13.54\) mg/h. Code that keeps adding \(R_0\) after \(T\) will fail this phase-specific check.

First-Order Absorption: Distinguish Input From Central Amount

For an oral dose represented by a depot and a central compartment, assume a first-order transfer from the depot to the central compartment and first-order elimination from the central compartment. Let \(A_g\) be depot amount and \(A_c\) central amount. For a dose \(D\), take \(A_g(0)=D\) and \(A_c(0)=0\). The model has no additional input after the initial depot amount is set.

Model: The equations are \(\frac{dA_g}{dt}=-k_aA_g\) and \(\frac{dA_c}{dt}=k_aA_g-k_eA_c\). The absorption flow into the central compartment is \(k_aA_g\); it is not the central amount \(A_c\).

Solving the depot equation gives \(A_g(t)=De^{-k_at}\). Substituting it into the central equation gives a linear differential equation with forcing term \(k_aDe^{-k_at}\). Multiplying by \(e^{k_et}\), integrating from zero to \(t\), and applying \(A_c(0)=0\) yields, when \(k_a\ne k_e\):

\frac{d}{dt}(e^{k_et}A_c) =k_aD e^{(k_e-k_a)t}, \qquad e^{k_et}A_c(t) =\frac{k_aD}{k_e-k_a}\left(e^{(k_e-k_a)t}-1\right), $$ $$ A_c(t)=\frac{Dk_a}{k_a-k_e} \left(e^{-k_et}-e^{-k_at}\right).

The initial conditions provide a particularly useful code check. At time zero, \(A_c(0)=0\), but the input flow is \(k_aA_g(0)=k_aD\), which is positive for positive \(k_a\) and \(D\). Consequently, the initial central derivative is \(dA_c/dt|_{t=0}=k_aD-k_e(0)=k_aD\). Confusing the zero central amount with zero central input reverses this interpretation.

Worked Example: Checking a Depot-to-Central Model

Consider a hypothetical oral model with \(D=100\) mg, \(k_a=1.20\) h\(^{-1}\), \(k_e=0.200\) h\(^{-1}\), and \(V=20.0\) L. The rates are unequal, so the derived expression applies directly. The initial states are \(A_g(0)=100\) mg and \(A_c(0)=0\) mg.

At \(t=1.00\) h, the expected central amount is:

A_c(1)=\frac{(100)(1.20)}{1.20-0.200} (e^{-0.200}-e^{-1.20}) =120(0.8187307531\ldots-0.3011942119\ldots) =62.10438494\ldots\ \text{mg}.

Therefore \(C(1)=62.10438494\ldots/20.0=3.105219247\ldots\) mg/L, or 3.105 mg/L. At time zero, the central amount is zero, while the central input is \(k_aA_g(0)=(1.20)(100)=120\) mg/h. The central amount begins increasing at 120 mg/h; that does not mean it jumps to 120 mg at time zero.

An implementation should therefore express the central derivative as absorption flow minus elimination flow, \(k_aA_g-k_eA_c\). A code expression that uses \(k_aA_c\) as the input instead would give zero input at time zero and would not match the model. Comparing both the initial derivative and a later concentration helps expose this kind of mismatch.

Diagnostics, Alternatives, and Limitations

A structural equation check is strongest when it uses independent features of the model rather than just one reported prediction. Useful checks include whether initial states match the stated conditions, whether derivatives have the expected signs, whether amounts remain nonnegative under the stated positive-input assumptions, whether a rate input switches off when intended, and whether the long-time behavior follows from the equations. For example, with no continuing input and a positive elimination rate, the one-compartment amount tends to zero. With continuing constant infusion and first-order elimination, the infusion-phase solution tends toward \(R_0/k\). These are checks of the stated equations, not universal facts about all pharmacokinetic models.

The main alternatives answer related but different questions. Comparing code with a written specification checks consistency with that specification; independently deriving equations checks whether the mathematical consequences of the specified structure are present. Running the code and comparing outputs with an independent analytical solution checks numerical behavior at selected conditions, but it may not reveal why a discrepancy occurs. Numerical simulation is useful when an analytical solution is unavailable, but then the reference must itself be independently constructed and checked. None of these approaches makes the others unnecessary.

Do not require textual identity between code and a derived formula. For instance, the infusion amount can be written as \(R_0(1-e^{-kt})/k\) or as \((R_0/k)(1-e^{-kt})\); these are algebraically equivalent when \(k\ne 0\). Conversely, two expressions that happen to agree at one time need not define the same trajectory. Compare equations, conditions, and multiple consequences. If the code uses a numerical solver, compare within a justified numerical tolerance and verify that the comparison uses the same inputs and evaluation times.

Common pitfall: “The central amount is zero at time zero, so the central input must also be zero.” This is wrong for the first-order absorption model. The amount is \(A_c(0)=0\), but the input is \(k_aA_g(0)=k_aD>0\) for positive \(k_a\) and \(D\). The correction is to distinguish a flow into a compartment from the amount already in that compartment.

Analytical checks have limits. They rely on the assumptions used in the derivation: for example, constant rate parameters, the specified compartment connections, and correctly stated initial conditions. A mismatch might reflect an implementation error, but it could also mean the reviewer compared different conditions or misunderstood a code construct. Record the equation, code location, parameters, initial states, and comparison time used so that another reviewer can reproduce the check. A matching result supports consistency for the conditions tested; it does not prove that every branch, event, or parameter value is correct.

Check Your Understanding

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

For \(dA/dt=-kA\) with \(A(0)=D\), what initial slope should the code imply?
In the infusion example, why is the post-infusion equation based on \(A(T)\) rather than resetting the amount to zero at \(T\)?
For the depot model, state the central amount and central input at time zero.
Why is agreement at a single time point insufficient to establish that two structural equations are equivalent?
Name one assumption that must be aligned when comparing an analytical equation with model code.

Solutions. (1) The initial slope is \(-kA(0)=-kD\). (2) Stopping the input does not remove the amount already in the compartment; the subsequent solution starts from the amount present at \(T\), namely \(A(T)\). (3) \(A_c(0)=0\), while the input is \(k_aA_g(0)=k_aD\), positive for positive dose and absorption rate. (4) Distinct functions can coincide at one time, so compare the equations or additional values and properties such as initial slope and limiting behavior. (5) Examples include matching the initial conditions, rate values, input-on/off phase, or evaluation time.