Adjoint Sensitivities for Many Parameters

Why Use an Adjoint?

A PBPK model can contain many parameters but be used to answer one particular quantitative question, such as the predicted concentration at a specified time or the total exposure over a dosing interval. Forward sensitivity equations calculate derivatives with respect to parameters by augmenting the state system with one sensitivity system per parameter. When the number of parameters is large and the quantity of interest is a single scalar, that approach can require many additional equations.

The adjoint method reverses the calculation. It first solves the PBPK state equations forward in time, then solves one auxiliary system backward in time. The resulting gradient contains derivatives with respect to all parameters. This is especially useful for local sensitivity calculations of a scalar endpoint; it does not make every sensitivity problem cheaper, because the cost depends on how many distinct outputs are needed.

Here \(\phi\) describes a terminal contribution, such as the amount in a compartment at time \(T\), and \(L\) describes a running contribution, such as the concentration accumulated over time. A terminal-only endpoint has \(L=0\). The partial derivatives \(f_x\) and \(f_\theta\) are Jacobian matrices, while \(\phi_x\), \(L_x\), and \(\lambda\) are treated as column vectors. The formulas below assume that \(f\), \(\phi\), and \(L\) are differentiable in the state and parameters, and that the state solution is sufficiently regular on the interval. Dosing events or other discontinuities require appropriate jump conditions, discussed below.

Adjoint Equation and Gradient

For the endpoint defined above, the adjoint is determined by integrating backward from \(T\):

Once the state trajectory and adjoint are available, the parameter gradient is

The term \(x_{0,\theta}\) is the Jacobian of the initial state with respect to the parameters. If the initial state is fixed independently of \(\theta\), this term is zero. The direct terms \(\phi_\theta\) and \(L_\theta\) matter when the endpoint depends explicitly on a parameter, even after the state trajectory is held fixed.

Why the Formula Holds

The derivation shows how the adjoint removes the need to calculate the full state-sensitivity matrix. Write \(\delta x(t)\) for a first-order state change caused by a parameter change \(\delta\theta\). Differentiating the state equation gives \(\delta\dot{x}=f_x\delta x+f_\theta\delta\theta\), and the endpoint change is \(\delta J=\phi_x^T\delta x(T)+\phi_\theta^T\delta\theta+\int_0^T(L_x^T\delta x+L_\theta^T\delta\theta)\,dt\). To eliminate the unknown \(\delta x(t)\) from the integral, add the zero quantity \(\int_0^T\lambda^T(f_x\delta x+f_\theta\delta\theta-\delta\dot{x})\,dt\). Integration by parts gives \(-\int_0^T\lambda^T\delta\dot{x}\,dt=-\lambda(T)^T\delta x(T)+\lambda(0)^T\delta x(0)+\int_0^T\dot{\lambda}^T\delta x\,dt\).

Choose \(\lambda(T)=\phi_x\) so the terminal-state variation cancels, and choose \(\dot{\lambda}=-f_x^T\lambda-L_x\) so the remaining terms involving \(\delta x(t)\) cancel at every time. Since \(\delta x(0)=x_{0,\theta}\delta\theta\), the terms left over are exactly the stated gradient formula. This argument requires differentiable state and endpoint functions; it is not a claim that the numerical solution has zero error.

Worked Calculations

Computational Cost, Events, and Checks

The main computational comparison is between the number of parameters \(p\) and the number of outputs \(q\). A forward sensitivity method propagates derivatives of the state with respect to all parameters, usually adding \(m\times p\) sensitivity values to the system. An adjoint method for one scalar output adds a backward system of dimension \(m\), plus the work of accumulating the gradient. This can be attractive when \(p\) is large and \(q\) is small. For \(q\) distinct scalar outputs, however, each generally has its own terminal condition and adjoint calculation. Forward sensitivities may be more efficient when there are few parameters and many outputs. Finite differences are often easy to implement but require repeated model runs, and their accuracy depends on the perturbation size.

The cost is not determined by parameter count alone. The state trajectory must be available during the backward solve, which may require storing it or reconstructing it with checkpointing. Evaluating \(f_\theta^T\lambda\) can itself be costly for a model with many parameters. The adjoint also requires a sufficiently accurate backward integration and a correct derivative of the implemented model; it does not repair an inaccurate state solution.

PBPK models often include scheduled bolus doses or other known state resets. At a reset \(x^+=R(x^-,\theta)\), the state is discontinuous and the adjoint must not simply be integrated through the event as if the state were smooth. For a smooth reset map, the adjoint satisfies \(\lambda^- = R_x^T\lambda^+\), and the parameter gradient receives the event contribution \(R_\theta^T\lambda^+\). Continuous infusion start or stop times may change the right-hand side and require integration to be split at the event time. More complicated event rules need their corresponding derivative conditions; ignoring an event can produce a plausible-looking but incorrect gradient.

Useful diagnostics include checking that the baseline state and endpoint agree with the ordinary PBPK solve, tightening solver tolerances to assess gradient stability, and comparing selected adjoint derivatives with central finite differences. For a parameter \(\theta_i\), the latter estimate is \([J(\theta_i+h)-J(\theta_i-h)]/(2h)\), with other parameters held fixed. Agreement should be assessed over reasonable perturbation sizes and solver tolerances; a single step size is not proof. Also check derivative units and signs, the treatment of initial-condition dependence, and whether all dosing events and explicit parameter terms have been included.

Common Pitfalls and Limits

• Using the wrong endpoint. A terminal concentration and an AUC have different endpoint definitions. The running term \(L\) changes the adjoint equation and must be included when the output is an integral.
• Assuming every parameter affects only the ODE. If \(\phi\), \(L\), or \(x_0\) depends directly on \(\theta\), its gradient contribution must be retained.
• Assuming one adjoint gives every output gradient. One adjoint corresponds to one scalar objective. Several outputs generally require several backward solves, though a weighted combination of outputs can be treated as one scalar objective when that is the quantity actually wanted.
• Ignoring numerical and event handling. State interpolation error, loose tolerances, missed resets, or incorrect event derivatives can undermine the result. Compare selected derivatives with an independent calculation and refine the numerical solution.

Adjoints calculate local derivatives at the specified parameter vector and model scenario. They do not by themselves establish that a parameter is identifiable, that a local effect persists over a broad parameter range, or that the model is adequate for a clinical prediction. Their role is narrower: efficiently differentiate a defined scalar model output, subject to the model equations, endpoint definition, and numerical implementation.