Debugging a Coupled QSP Model Step by Step

Debugging the Model, Not Just the Solver

A coupled QSP simulation can fail in a way that looks biologically plausible: drug concentration falls, a biomarker changes, and the solver returns a smooth curve. Smoothness is not evidence that the model is correct. A reliable debugging sequence reproduces the run, checks the model’s states and inputs, and finds the first place where the computed behavior conflicts with the stated equations.

We will use a small drug-exposure and biomarker-turnover model. Its size makes every term inspectable, while its coupling captures a common QSP structure: drug exposure changes a process that controls biomarker production. The goal is not to tune a solver until the plot looks reasonable; it is to identify and test the equation or input responsible for the failure.

Model and assumptions: Let \(C(t)\) be drug concentration in mg/L and \(B(t)\) be biomarker concentration in biomarker units/L. Time \(t\) is in hours. After a single bolus at \(t=0\), with no further doses or other drug input, the model is

\frac{dC}{dt}=-k_e C,\qquad \frac{dB}{dt}=k_{\mathrm{out}}\left[B_0\left(1-I(C)\right)-B\right], \qquad I(C)=\frac{I_{\max}C}{IC_{50}+C}.

Here \(k_e\) and \(k_{\mathrm{out}}\) have units h\(^{-1}\); \(B_0\) and \(B\) have the same biomarker units; \(I_{\max}\) is dimensionless; and \(IC_{50}\) has units mg/L. Assume nonnegative concentrations, positive rate constants and \(IC_{50}\), \(0\leq I_{\max}\leq1\), and a drug effect that reduces biomarker production. The bolus is represented by \(C(0)=D/V\), where dose \(D\) and volume \(V\) use compatible units. The baseline condition is \(B(0)=B_0\).

The state vector is ordered as \((C,B)\), so the right-hand side must return the two derivatives in that same order. A useful first check is units: \(dC/dt\) must be in mg/L/h, and \(dB/dt\) must be in biomarker units/L/h. The inhibition \(I(C)\) is dimensionless, so \(B_0(1-I(C))\) has the units of \(B\), as required.

Derive Checks Before Running the Full Model

The equation for \(C\) has an analytic solution. Separating variables and applying \(C(0)=C_0\) gives \(C(t)=C_0e^{-k_et}\). This gives a reference for exposure independent of the biomarker equation. At any fixed concentration \(C\), the biomarker target is \(B^*(C)=B_0(1-I(C))\). Subtracting this target from \(B\) gives

\frac{d}{dt}\left(B-B^*(C)\right)=-k_{\mathrm{out}}\left(B-B^*(C)\right)

when \(C\) is held fixed. Thus the fixed-exposure equilibrium is \(B^*(C)\), and deviations from it decay at rate \(k_{\mathrm{out}}\). In the full model, \(C\) changes with time, so this fixed-exposure result is a limiting or unit-test check, not a claim that \(B\) instantly equals its moving target.

At the start of the simulation, calculate the derivatives directly from the equations. These initial slopes are often more diagnostic than a long trajectory: a wrong sign, swapped state, or incorrect initial concentration can be visible immediately. For the values below, the drug must initially decline and the biomarker must initially decline because production is inhibited while \(B(0)=B_0\).

Worked Example: Reproducing a Sign Error

Suppose \(k_e=0.5\) h\(^{-1}\), \(C(0)=4\) mg/L, \(k_{\mathrm{out}}=0.2\) h\(^{-1}\), \(B_0=B(0)=100\) biomarker units/L, \(I_{\max}=0.8\), and \(IC_{50}=2\) mg/L. A failed implementation uses \(dB/dt=k_{\mathrm{out}}[B-B_0(1-I(C))]\), reversing the sign of the biomarker turnover term. Reproduce the failure at the initial time.

First evaluate the inhibition and target:

I(4)=\frac{0.8(4)}{2+4}=0.5333,\qquad B^*(4)=100(1-0.5333)=46.67\ \text{biomarker units/L}.

The drug slope is \(dC/dt=-0.5(4)=-2.000\) mg/L/h. The correct biomarker slope is \(0.2(46.67-100)=-10.67\) biomarker units/L/h. The faulty equation instead gives \(0.2(100-46.67)=+10.67\) biomarker units/L/h: it predicts that the biomarker initially rises despite the model’s assumed inhibitory effect.

For a one-step reproduction, use explicit Euler with \(\Delta t=0.1\) h. The correct update is \(B_1=100+0.1(-10.67)=98.93\) biomarker units/L. The faulty update is \(B_1=100+0.1(10.67)=101.1\) biomarker units/L. This is a targeted reproduction: the divergence occurs in the first state update, so changing the solver or its tolerances would not repair the equation.

Isolate the Fault in a Reproducible Sequence

Keep the failed run’s parameter values, state ordering, initial conditions, and dosing inputs fixed while debugging. Otherwise, changing several things at once can hide which change mattered. Work from the simplest checks toward the full simulation: inspect inputs and units, evaluate the right-hand side at a known state, then compare a short run with an independent reference.

Reproduce and record the failure.
Run the same initial state, parameters, and input schedule. Record the first warning or the first time and state at which the output becomes implausible. Save the state order and units alongside the numbers.

Check the inputs and initial state.
Verify that the bolus has been converted to concentration, that the biomarker starts at its stated baseline, and that all parameters are in the units expected by the equations.

Evaluate each derivative separately.
At a hand-checkable state, calculate the exposure and biomarker derivatives from the equations. Compare each result with the derivative returned by the implemented right-hand side, in the documented state order.

Repair one identified cause and retest.
Change only the faulty equation or input. Repeat the same small test, then run targeted boundary and trajectory checks before comparing the repaired full run with the original failure.

A model-level test should be independent enough to catch the suspected error. Recomputing the same faulty expression in a test does not provide an independent check. For example, calculate the expected initial slopes by hand, or compare the exposure state with its analytic solution. Earlier material on translating ODE systems into right-hand-side functions and mapping biological states to indexed vectors is useful here: the implemented output must follow the model’s stated state order.

Worked Example: Separating a Dose-Conversion Error From an ODE Error

Use the same parameters as above, but let the administered dose be \(D=200\) mg and the apparent distribution volume be \(V=50\) L. The correct initial concentration is \(C(0)=D/V=4\) mg/L. Suppose the implementation mistakenly assigns the dose value, 200, directly to the concentration state.

At the correct initial concentration, \(I(4)=0.5333\) and \(dB/dt=-10.67\) biomarker units/L/h, as calculated above. With the faulty state value \(C(0)=200\) mg/L, the implemented model instead evaluates

I(200)=\frac{0.8(200)}{2+200}=0.7921,\qquad B^*(200)=100(1-0.7921)=20.79\ \text{biomarker units/L},

so \(dB/dt=0.2(20.79-100)=-15.84\) biomarker units/L/h. Both runs have a negative initial biomarker slope, so a direction-only check would miss the problem. The magnitude and exposure state reveal it. Inspecting the initial state and checking \(D/V\) isolates the input conversion; the turnover equation itself is unchanged. In this example, the concentration error also makes the initial drug derivative \(-0.5(200)=-100\) mg/L/h rather than \(-2.000\) mg/L/h.

The corrective test is direct: assert that the bolus produces \(C(0)=200\text{ mg}/50\text{ L}=4\) mg/L, then re-evaluate both initial derivatives. A plausible biomarker curve is not a substitute for this input check, because the wrong concentration can still drive the model in the expected direction.

Confirm the Repair With Targeted Tests

After finding a likely fault, test the corrected model at states where the expected behavior is known. These tests check different parts of the implementation and should be interpreted together. A failed equilibrium test points toward the production or turnover expression; a failed analytic exposure comparison points toward the exposure equation, state order, or input; and a failed initial-slope check points toward the local right-hand side or its inputs.

Worked Example: Testing Baseline, Limiting Effect, and Exposure

For the parameter values in the first example, test three properties of the repaired implementation. First, when \(C=0\) and \(B=B_0=100\) biomarker units/L, inhibition is zero and the biomarker derivative must be zero:

I(0)=0,\qquad \frac{dB}{dt}=0.2[100(1-0)-100]=0\ \text{biomarker units/L/h}.

Second, at a very large concentration, \(I(C)\) approaches \(I_{\max}=0.8\). At \(B=100\), the limiting biomarker derivative is therefore \(0.2[100(1-0.8)-100]=-16.00\) biomarker units/L/h. This checks both the maximum-effect limit and the sign of the production effect.

Third, for \(C(0)=4\) mg/L, the analytic exposure at \(t=2\) h is \(C(2)=4e^{-0.5(2)}=1.472\) mg/L. A numerical run should approach this value as numerical error is reduced. The test applies specifically to \(C\); it does not assert that the biomarker follows the same exponential curve. Together, the baseline, high-exposure, and analytic-exposure checks probe distinct behaviors without requiring a full reference trajectory for \(B(t)\).

For the biomarker, also inspect whether the trajectory moves in a direction consistent with its current target. If \(B>B^*(C)\), then \(dB/dt<0\); if \(B<B^*(C)\), then \(dB/dt>0\). At fixed \(C\), the turnover equation pulls \(B\) toward \(B^*(C)\), rather than pushing it away. These directional checks help identify a sign error, but they do not validate every aspect of a trajectory: the dose conversion and parameter units still need separate checks.

Diagnostics, Alternatives, and Limitations

Warnings and implausible outputs are clues, not diagnoses. A solver failure may arise from a coding error, inconsistent units, inappropriate initial conditions, or numerical difficulty. Conversely, a solver can successfully integrate an incorrectly specified equation. Look for the first divergence and use the model structure to narrow the cause.

Wrong initial slope: evaluate the right-hand side at the initial state. Check signs, state order, and whether the computed derivative has the expected units.
Unexpected initial concentration: compare the implemented state with dose divided by volume, and confirm that the bolus is represented only once.
Correct start, implausible later behavior: inspect the changing input and each state’s derivative at the first time the trajectory departs from expectation. Do not infer the cause from the final value alone.
Negative or non-finite values: determine whether they are compatible with the model assumptions. In this model with nonnegative initial states and stated parameter constraints, negative concentrations warrant investigation; numerical output alone does not establish a biologically meaningful negative concentration.
Different results after a solver change: compare the same model and inputs under numerical convergence checks. A solver setting can affect numerical error, but it cannot correct a wrong sign or dose conversion.

There are several ways to investigate a failure. Looking at the final biomarker value is quick, but it can conceal compensating errors. Inspecting plots helps reveal timing and direction, but visual agreement is a weak test of equations. Reproducing individual derivatives, checking units and initial states, and using analytic or limiting-case references are more targeted. A full regression comparison is valuable after the repair, but it should not replace a test that identifies why the original result was wrong. Numerical convergence checks address sensitivity to numerical settings; they do not establish that the model equations or inputs are correct.

The model used here is deliberately simple. It assumes first-order drug elimination, a direct concentration-dependent effect on biomarker production, and a single bolus input. It does not describe distribution between compartments, indirect mechanisms beyond the stated turnover equation, parameter uncertainty, or additional dosing. A successful debugging test establishes that specific computations behave as expected under these assumptions; it does not establish that the assumptions adequately describe a particular biological system.

Key takeaway: Debug a coupled model by finding its first incorrect state or derivative, not by adjusting settings until the output looks acceptable. Reproduce the run, test inputs and equations separately, repair one identified cause, and verify the repair with independent checks.

Check Your Understanding

Use the model and assumptions above. Each problem includes a solution so you can check the reasoning.

If \(C=0\) and \(B=80\) biomarker units/L, with \(B_0=100\) and \(k_{\mathrm{out}}=0.2\) h\(^{-1}\), what is \(dB/dt\), and which direction should \(B\) initially move?
A simulation uses \(D=120\) mg and \(V=30\) L but sets \(C(0)=120\). What initial concentration should be used, and which test most directly identifies the error?
For \(C=2\) mg/L, \(I_{\max}=0.8\), and \(IC_{50}=2\) mg/L, calculate \(I(C)\) and \(B^*(C)\) when \(B_0=100\) biomarker units/L.
Why does matching the analytic \(C(t)\) solution not, by itself, verify the biomarker equation?
A biomarker rises initially when \(C>0\) and \(B=B_0\). Name a direct model check before changing the solver.

Solutions. (1) With no drug, \(I(0)=0\), so \(dB/dt=0.2(100-80)=+4.000\) biomarker units/L/h; \(B\) initially rises toward baseline. (2) The correct value is \(120/30=4\) mg/L. Check the dose-to-concentration conversion and assert the initial exposure state before integrating. (3) \(I(2)=0.8(2)/(2+2)=0.4000\), so \(B^*(2)=100(1-0.4)=60.00\) biomarker units/L. (4) The analytic solution tests the exposure equation and its implementation, but not the separate biomarker turnover expression or the coupling through \(I(C)\). (5) Evaluate \(I(C)\), the target \(B_0(1-I(C))\), and the implemented initial biomarker derivative; with the stated inhibitory model at baseline, that derivative should be negative.