Second-Order Local Sensitivity and Interactions

Why First-Order Rankings Can Miss Local Behavior

A first-order sensitivity describes how a PBPK output changes for a small change in one parameter while the other parameters are held fixed. It is a local slope. But a slope alone does not describe whether that slope itself changes nearby, or whether the effect of changing one parameter depends on the value of another. These features are local curvature and parameter interaction, respectively.

Second-order local sensitivity analysis adds the relevant second derivatives to the first-order picture. It can reveal that an output bends noticeably even when its slope at the reference parameter set is small, or that two parameters jointly affect an output in a way that cannot be represented by adding their separate first-order changes. The analysis remains local: it characterizes the model near a specified parameter set, not across the full plausible parameter space.

Parameter Coordinates and the Second-Order Approximation

PBPK parameters often have different units and scales, so derivatives with respect to raw parameters can be difficult to compare. For positive parameters, a useful local coordinate is the log-relative change \(u_i=\ln(\theta_i/\theta_{i0})\). The reference value corresponds to \(u_i=0\), and a small \(u_i\) represents an approximately proportional change in \(\theta_i\). For example, \(u_i=\ln(1.1)\) represents a 10% increase. Define \(g(\mathbf{u})=f(\boldsymbol{\theta})\) using this coordinate system. Derivatives below refer to these log-relative coordinates unless stated otherwise.

The gradient contains the first derivatives, and the Hessian contains the second derivatives at the reference point. If \(g\) is twice continuously differentiable near that point, its local second-order approximation is:

The approximation error is smaller than second order as \(\mathbf{u}\) approaches zero, under the stated smoothness condition. For one parameter, the second-order contribution is \(\frac{1}{2}H_{ii}u_i^2\). For two distinct parameters \(i\) and \(j\), the two corresponding terms in the double sum combine to \(H_{ij}u_i u_j\) when the mixed partial derivatives are equal. Thus the mixed derivative quantifies how the joint change differs from the sum of the separate local slopes, to second order.

A short derivation follows by considering a straight path from the reference point to the perturbed point: define \(\phi(t)=g(t\mathbf{u})\) for \(0\leq t\leq1\). The chain rule gives \(\phi'(0)=\sum_i (\partial g/\partial u_i)u_i\), and differentiating again gives \(\phi''(0)=\sum_i\sum_j (\partial^2g/\partial u_i\partial u_j)u_i u_j\). The one-dimensional second-order Taylor expansion of \(\phi(1)\) about zero is therefore exactly the displayed approximation, with a remainder that is \(o(\|\mathbf{u}\|^2)\) when the second derivatives are continuous at the reference point.

Worked Uses of Curvature and Mixed Partials

Estimating Mixed Derivatives from Model Runs

When analytic derivatives are unavailable, a mixed partial can be estimated by evaluating the model at four nearby points. With step sizes \(h\) and \(k\) in two log-relative coordinates, the centered finite-difference estimate is:

Here all coordinates other than \(i\) and \(j\) are held at their reference values. This formula uses changes in both directions so that constant and first-order terms cancel. It is an approximation for a smooth model; its accuracy depends on the step sizes and numerical error in the model evaluations.

Assumptions, Diagnostics, and Failure Modes

The Taylor approximation and symmetric mixed-partial interpretation require a sufficiently smooth output in the neighborhood being examined. If the output has a discontinuity, a sharp threshold, or a change in its definition, local derivatives may not exist or may be a poor summary. PBPK models may include conditional rules, dose or event handling, or numerical solver behavior that creates apparent irregularities. A derivative calculated across such a feature should not automatically be interpreted as biological curvature.

Numerical estimates should be checked by repeating them with more than one reasonable step size. If a finite-difference estimate changes greatly as the step is reduced, possible causes include truncation error at large steps, roundoff or solver error at very small steps, or a genuinely narrow region of changing curvature. Check that the four model runs use the same model settings, output definition, and solver tolerances. Where possible, compare finite differences with derivatives from an analytic or sensitivity-equation implementation.

How This Complements Other Sensitivity Approaches

First-order local sensitivity rankings summarize slopes at a reference parameter set. They are simpler to calculate and interpret when the aim is a small, one-parameter perturbation. Adding second derivatives is useful when curvature or joint local change matters, but it requires more derivatives or model runs and does not turn a local analysis into a global one.

For PBPK ODE models, forward sensitivity equations can be extended to compute second-order sensitivities, while finite differences can estimate them from repeated model evaluations. The former can be efficient when derivatives are implemented reliably; the latter is straightforward to apply but can be sensitive to step size and solver error. Neither approach removes the need to define the output, reference set, and parameter coordinates clearly.

The Hessian is also coordinate-dependent. A mixed partial computed using raw clearance and raw volume parameters is not numerically interchangeable with one computed in log-relative coordinates. Even its sign can change under nonlinear reparameterization in some circumstances. Further, a nonzero mixed partial is evidence of local mathematical dependence in the model response; by itself it does not establish a biological mechanism, parameter correlation, or causal interaction.

Common Mistake and Practical Limitations

• Mistake: “The mixed derivative is positive, so increasing both parameters must increase the output.” This ignores the signs and sizes of both perturbations, the first-order terms, and diagonal curvature. The mixed contribution is \(H_{ij}u_i u_j\), so its sign depends on both \(H_{ij}\) and the direction of change.
• Correction: Evaluate the full second-order expression for the perturbation of interest, then verify it against model runs when the change is not very small.
• Limitation: A Hessian describes local behavior at one reference set. It may not remain informative after large changes, in another physiological scenario, or for a different endpoint.

For reporting, include the output and units, reference parameter set, coordinate system, derivative units, and whether values were analytic or numerically estimated. If finite differences were used, report step sizes and a stability check. These details are necessary to distinguish a reproducible local result from a number that depends on an undocumented computational choice.