Stable Estimation of Stiff TMDD ODE Models

Why Numerical Stability Matters in TMDD

Target-mediated drug disposition (TMDD) models can combine rapid association and dissociation with much slower target turnover, drug elimination, or dosing intervals. A numerical solver must resolve the rapid changes without taking steps so large that it misses them. At the same time, it must continue to integrate accurately across the slower parts of the profile. This combination can make the system stiff: stable and accurate computation may require many small steps from an explicit solver, even when the quantities of scientific interest change slowly.

This is not merely a computational inconvenience. Parameter estimation repeatedly solves the differential equations at trial parameter values. If the solver misses a fast transient, accumulates substantial integration error, or fails unpredictably, the resulting predictions—and therefore the likelihood being optimized—can be inaccurate or irregular. A fitted parameter can then reflect numerical behavior rather than the data.

For illustration, consider a concentration-based TMDD model with free drug \(D\), free target \(R\), and drug–target complex \(C\), all measured in nM. Assume the pools share a compartment and effective volume, input \(u(t)\) is expressed as a concentration per hour, and the rate constants and synthesis rate are constant between dosing events. A basic full model is

Here \(k_{\mathrm{on}}\) has units nM\(^{-1}\)h\(^{-1}\); \(k_{\mathrm{off}}, k_{\mathrm{el}}, k_{\mathrm{deg}}\), and \(k_{\mathrm{int}}\) have units h\(^{-1}\); and \(k_{\mathrm{syn}}\) has units nM/h. Binding removes free drug and free target while adding complex. Dissociation reverses that transfer; complex internalization removes complex. This model assumes nonnegative initial concentrations and rates, and it treats the stated processes as the mechanisms represented by these equations. More elaborate compartment structures may require additional states and equations.

Where the Fast Time Scale Comes From

A useful local explanation comes from considering binding and dissociation over a short interval, while treating the slower processes as negligible. If total drug and total target are approximately conserved during that interval, then \(D=D_{\mathrm{tot}}-C\) and \(R=R_{\mathrm{tot}}-C\). The complex obeys

At a local state with free concentrations \(D\) and \(R\), the derivative of the right-hand side with respect to \(C\) is \(-k_{\mathrm{on}}(D+R)-k_{\mathrm{off}}\). Thus a small displacement from the local binding balance decays at approximately this rate. When internalization also acts on the fast interval, a useful local approximation includes its removal rate:

This is a time-scale guide, not a complete eigenvalue analysis of the full model. The full Jacobian also includes turnover, elimination, input, and coupling to other states. Stiffness becomes plausible when the fast decay time is far shorter than important slow time scales, such as \(1/k_{\mathrm{deg}}\) or \(1/k_{\mathrm{el}}\). The separation can change over time because \(D\) and \(R\) change.

Solver Tolerances and State Scales

Adaptive ODE solvers estimate local numerical error and adjust step size. A common scale for the allowed local error in state \(y_i\) is

where \(\mathrm{rtol}\) is a dimensionless relative tolerance and \(\mathrm{atol}_i\) is an absolute tolerance for state \(i\), expressed in that state’s units. This is a practical error-control scale, not a guarantee that the global error in a predicted concentration is below the same threshold. The solver’s error estimate, method, and behavior across steps all matter.

A single absolute tolerance can be inappropriate when states have different typical magnitudes. If a target concentration is \(10^{-3}\) nM while drug concentrations are around \(1\) nM, an absolute tolerance adequate for drug may be too loose for target. Conversely, demanding extremely small absolute errors for every state can waste computation or expose round-off limitations. Set tolerances with the state units, expected magnitudes, and the precision needed for the model predictions in mind. Rescaling states or using state-specific absolute tolerances can help.

Relative tolerance alone is also insufficient near zero: when \(y_i\) is close to zero, \(\mathrm{rtol}|y_i|\) is close to zero, so the absolute tolerance controls the error scale. This matters for concentrations that become very small, especially if predictions or likelihoods depend on their accurate values.

Estimation Strategies That Preserve the Likelihood

For observations \(y_j\) with model predictions \(f_j(\theta)\), a likelihood calculation depends on the predictions at parameter vector \(\theta\). For example, under independent Gaussian errors with known standard deviations \(\sigma_j\), the part of \(-2\log L\) that depends on the predictions is

A small change in a prediction can matter if it is large relative to its observation error. Solver accuracy therefore needs to be judged in relation to the likelihood’s sensitivity, not just by whether plotted curves look similar. A solver failure at a trial parameter value is not a valid prediction. It should be handled explicitly by the estimation procedure, rather than silently omitting the affected observations or substituting an arbitrary likelihood contribution.

A practical estimation workflow uses a solver designed for stiff systems when the model exhibits fast and slow dynamics. Implicit methods such as backward differentiation formulas or Radau methods are standard choices; explicit methods can be effective when the system is not stiff or when the relevant interval is short, but may require prohibitively small steps in stiff regions. Supplying a Jacobian or its sparsity pattern, when supported and correctly specified, can improve efficiency for larger systems.

Parameter scaling matters as well. TMDD parameters can differ greatly in magnitude and units. Optimizers often behave better when positive parameters are represented on a log scale, such as optimizing \(\log k_{\mathrm{on}}\) rather than \(k_{\mathrm{on}}\), and when the parameterization avoids implausible or numerically problematic regions. This does not remove the need for scientifically appropriate bounds, initial values, or identifiability assessment. It helps the optimizer navigate scale; it cannot create information absent from the data.

Diagnostics, Alternatives, and Limitations

Warning signs include repeated solver step-size reductions, excessive step counts, failed steps, non-finite predictions, negative concentrations when the model should remain nonnegative, and large changes in predictions or estimates after modest tolerance tightening. These signs call for diagnosis rather than automatic changes to the scientific model. Check units and rate scales, input discontinuities such as dosing events, initial conditions, state scaling, and whether a solver is appropriate for stiffness. A very small tolerance is not by itself proof of accuracy: it may increase cost without fixing an incorrectly coded equation or an unsuitable parameterization.

The main numerical alternative to an implicit stiff solver is an explicit solver. Explicit methods are often simpler and can be fast for nonstiff systems, but stability restrictions may force tiny steps in the presence of fast binding. A second alternative is to use a reduced TMDD model or a rapid-binding approximation. That changes the model’s mathematical assumptions and may change its predictions; it is not simply a faster solver setting. Use such a reduction only when its assumptions are justified for the scientific question and data. The full-versus-reduced model choice is distinct from the numerical task of accurately solving whichever equations are chosen.

Numerical checks also have limits. Agreement between two solvers does not prove the model is correct, and stable optimization does not establish that parameters are identifiable or biologically interpretable. Conversely, an optimizer’s failure can arise from poor scaling or integration failures rather than from a scientifically implausible model. Numerical reliability is one necessary part of a defensible fit, alongside model assumptions, data information, and appropriate evaluation of fit.