Stable Numerical Estimation of Full TMDD ODEs

Why Numerical Stability Matters in Full TMDD Estimation

A full target-mediated drug disposition (TMDD) model can contain processes that evolve on very different time scales. Binding may be rapid, while target turnover and drug elimination are slower. An integration method that takes steps suitable for the slow dynamics can miss the rapid ones; a likelihood optimizer can then respond to numerical error as if it were information about the parameters. The result may be unstable estimates, unreliable gradients, or apparent optimizer convergence at a parameter value that changes when solver settings change.

The task here is not to choose between a full TMDD model and a reduced model. It is to make the numerical solution of the full system, and the likelihood calculations based on it, sufficiently accurate and reproducible. The essential checks are to control integration error, represent dose events correctly, and verify that parameter derivatives remain stable at the accuracy used for fitting.

Full TMDD Equations and the Source of Stiffness

For clarity, write \(C(t)\) for free drug concentration, \(R(t)\) for free target concentration, and \(B(t)\) for drug–target complex concentration. Assume these states use compatible concentration units. Let \(u(t)\) be drug input expressed as amount per time, \(V\) the volume, \(k_e=CL/V\) the first-order elimination rate for free drug, \(k_{\mathrm{on}}\) and \(k_{\mathrm{off}}\) the binding and dissociation rate constants, \(k_{\mathrm{syn}}\) and \(k_{\mathrm{deg}}\) the target synthesis and degradation rates, and \(k_{\mathrm{int}}\) the complex internalization rate. Under mass-action kinetics, the equations are

These equations presume the stated reaction scheme and parameterization, nonnegative initial states, and a correctly specified input function. In the absence of drug and complex, the target equilibrium is \(R_0=k_{\mathrm{syn}}/k_{\mathrm{deg}}\). The dynamics can be stiff when, for example, binding and dissociation rates are much faster than target turnover. Stiffness is a numerical property: it means that stable integration may require small steps even where the solution itself changes slowly. It does not, by itself, mean the model is wrong.

An explicit solver is often effective for nonstiff systems, but can require prohibitively small steps on a stiff one. Implicit solvers, such as BDF or Radau methods, are common choices for stiff systems. Solver choice is not a substitute for checking accuracy: even a suitable stiff solver can return inadequate results if tolerances are too loose, state scales are poorly handled, or dosing events are represented incorrectly.

Tolerances, Sensitivities, and the Likelihood

An adaptive solver estimates local integration error and adjusts its step size to keep that error within a tolerance. A common control uses a scale such as \(\mathrm{atol}_i+\mathrm{rtol}|x_i|\) for state \(x_i\): \(\mathrm{rtol}\) sets relative accuracy, while \(\mathrm{atol}_i\) provides an absolute scale when a state is near zero. In a TMDD system, concentrations of \(C\), \(R\), and \(B\) may differ substantially, so a single absolute tolerance may be inappropriate. Use state-specific absolute tolerances when needed, with units and magnitudes chosen to resolve scientifically relevant changes in each state.

The solver's local error control does not guarantee that the final likelihood or its gradient is accurate enough for estimation. Small integration errors may accumulate or affect predictions at observation times. Therefore, check stability of predicted concentrations and the objective function as tolerances are tightened. The appropriate tolerance is one for which the quantities used in fitting are stable to a level that is negligible relative to the data's measurement error and the precision required for the estimates.

For a differentiable ODE \(dx/dt=f(x,q,t)\), differentiating both sides with respect to \(q_j\) gives the sensitivity equation

This follows from the chain rule: the parameter affects the derivative both through the current state \(x\) and directly through \(f\). For the full TMDD equations, the state Jacobian is

For example, if \(q_j=\log(k_{\mathrm{on}})\), the direct derivative \(\partial f/\partial q_j\) is \((-k_{\mathrm{on}}CR,-k_{\mathrm{on}}CR,k_{\mathrm{on}}CR)^\mathsf{T}\). Sensitivities also need initial values: if an initial state depends on a parameter, differentiate that initial condition as well. For instance, if \(R(0)=k_{\mathrm{syn}}/k_{\mathrm{deg}}\), its sensitivity to either target-turnover parameter is generally not zero.

At each observation time, state sensitivities pass through the observation model by the chain rule. They can then contribute to derivatives of the likelihood. Sensitivity equations are one way to obtain gradients; finite differences estimate a gradient by perturbing parameters and rerunning the ODE solver. Automatic differentiation is another option when it correctly differentiates the numerical solution and event handling. Each approach still depends on accurate state solutions and correct treatment of inputs.

Discontinuous Dosing Inputs and Event Handling

Dose records are inputs to the ODE, not ordinary smooth functions. A bolus dose at time \(t_d\) changes the drug state instantaneously; it does not cause an instantaneous change in free target or complex under the model above. For a dose amount \(D\), the drug concentration jump is

A constant-rate infusion that starts or stops changes \(u(t)\) discontinuously, but the states remain continuous at the switch unless the record also includes a bolus. Integrate up to each event, apply the appropriate jump if there is one, and restart integration with the new input. Allowing a solver step to span an unrepresented input change can smear the dose across time and corrupt predictions and gradients.

Gradient Checks and Common Failure Modes

A likelihood optimizer uses changes in the objective to decide how to update parameters. If a gradient is noisy because repeated ODE solutions have inconsistent numerical error, an optimizer may take unproductive steps, stop prematurely, or report convergence for the wrong reason. Parameter transformations help enforce constraints—for example, optimizing \(q=\log(k_{\mathrm{on}})\) keeps \(k_{\mathrm{on}}>0\)—but do not remove numerical error. Gradients should be checked at plausible parameter values and, when possible, near the reported solution.

Common warning signs include predictions that change materially under tighter tolerances, a large number of rejected solver steps, negative concentrations beyond negligible rounding error, abrupt jumps at times without a dose, and gradients that change substantially with small changes in finite-difference step size. These signs call for investigation, not automatic parameter adjustment. Check state scales and absolute tolerances, event times and units, solver suitability, parameter transformations, and the derivative calculation. A negative state can also reflect an unsuitable model or parameter region; clipping it to zero without justification changes the equations and can introduce a nondifferentiable operation.

A practical comparison is to fit or evaluate the same objective with a stiff solver at progressively tighter tolerances, then compare the predictions, objective values, and gradients. A nonstiff explicit method may provide a useful independent check if it can solve the problem accurately, but failure or excessive computation with that method is not evidence that the model is invalid. Likewise, finite-difference gradients are useful checks but can be sensitive to both perturbation size and solver error. Sensitivity equations are often efficient when many likelihood evaluations are needed, but require correct derivative equations and initial sensitivities.

Numerical stability is necessary but not sufficient for sound inference. Stable integration cannot resolve structural non-identifiability, inadequate data, or misspecification of the observation model. Conversely, unstable estimates that move as tolerances change should not be interpreted biologically until numerical accuracy is addressed. The appropriate standard is not a particular solver setting in isolation, but reproducibility of the quantities used for estimation at a precision adequate for the analysis.