Why Differentiate the PBPK Equations Directly?
A PBPK model often represents organ amounts, tissue concentrations, or both as a system of differential equations. A local sensitivity asks how a model state changes when a parameter changes slightly around a specified value. Finite-difference calculations approximate this change by running the model at nearby parameter values. Forward sensitivity equations instead differentiate the ODE system itself and integrate the resulting derivatives alongside the states.
This approach is especially useful when a smooth model must be evaluated repeatedly at specific parameter values, or when numerical differences are unstable because small perturbations are obscured by solver error. It requires derivatives of the model equations and adds one sensitivity state for each state-parameter pair. The result is a local derivative, not a global measure of parameter influence.
Deriving the Forward Sensitivity Equation
Consider a smooth state equation \(dx/dt=f(t,x,\theta)\), with initial state \(x(0,\theta)=x_0(\theta)\). Assume that \(f\) is differentiable in the state and parameter over the interval of interest, and that the solution remains well defined there. The parameter is held fixed in time. Differentiating the state equation with respect to \(\theta\), and applying the chain rule to the right-hand side, gives the forward sensitivity equation.
Here, \(\partial f/\partial x\) is the Jacobian of the state derivatives with respect to the states, and \(\partial f/\partial\theta\) is the direct parameter derivative while holding the states fixed. The first term propagates sensitivity already present in the state vector; the second adds the parameter's direct effect on the rates of change. If \(\theta\) is one of several parameters, the equation is solved for each parameter separately, using the corresponding column of derivatives.
The initial condition matters. A parameter-dependent initial amount has a nonzero initial sensitivity even if the parameter does not appear in the subsequent rate equations. Conversely, an initial state independent of \(\theta\) has \(S(0)=0\), but the sensitivity can become nonzero later through \(\partial f/\partial\theta\).
Continuous Inputs and Dosing Events
A continuous input belongs in the state equation. If an amount \(A\) follows \(dA/dt=-kA+R(\theta)\), then its sensitivity equation includes the derivative of the input rate, \(dR/d\theta\), as well as the state-propagation term. Omitting the input derivative would treat the infusion rate as fixed even when it is the parameter being varied.
A bolus is different: it changes the state instantaneously rather than acting as a smooth rate over time. Write the event as \(x^+=g(x^-,\theta)\), where \(x^-\) and \(x^+\) are the states immediately before and after the dose. At a fixed, parameter-independent dosing time, differentiating this jump gives the sensitivity update.
For a dose added to one amount, \(A^+=A^-+F D\), the jump derivative with respect to bioavailability \(F\) is \(S_A^+=S_A^-+D\). A simulation should integrate up to the dose time, apply both state and sensitivity jumps, and continue from the post-dose values. A smooth ODE solver should not be asked to smear an instantaneous dose across a time step.
If the dosing time itself depends on the parameter, the event timing also contributes. For event time \(\tau(\theta)\), the fixed-clock-time sensitivity just after the event is \(S^+=g_xS^-+g_\theta+(g_xf^- - f^+)\,d\tau/d\theta\). The extra term accounts for the time spent under the pre-event and post-event dynamics changing as the event shifts. This correction is unnecessary when dosing times are fixed.
Worked Examples
Worked Example: Sensitivity to a Constant Infusion Rate
Suppose a hypothetical one-compartment amount \(A(t)\), in mg, obeys \(dA/dt=-kA+R\), with \(A(0)=0\). Let \(k=0.20\ \mathrm{h}^{-1}\) and \(R=4.00\ \mathrm{mg/h}\). Find the sensitivity to \(R\) and use it to calculate the amount at \(t=2.00\) h.
Derive the equation. Here \(f=-kA+R\), so \(\partial f/\partial A=-k\) and \(\partial f/\partial R=1\). Thus \(dS_R/dt=-kS_R+1\), with \(S_R(0)=0\), because the initial amount does not depend on \(R\). Its solution is \(S_R(t)=(1-e^{-kt})/k\).
Substitute and calculate. At 2.00 h, \(S_R(2)=(1-e^{-0.20(2.00)})/0.20=(1-e^{-0.4})/0.20=1.6484\ \mathrm{h}\). The amount is \(A(2)=R S_R(2)=4.00(1-e^{-0.4})/0.20=6.5936\ \mathrm{mg}\), or approximately \(6.594\ \mathrm{mg}\).
The sensitivity has units of hours: a local increase of \(1\ \mathrm{mg/h}\) in the infusion rate corresponds to an increase of about \(1.6484\ \mathrm{mg}\) in the amount at 2.00 h. This derivative is exact for the stated linear model; in a nonlinear PBPK model, the local derivative need not remain constant over a large parameter change.
Worked Example: Bioavailability in a Parameter-Dependent Bolus
Consider a hypothetical amount \(A\) that receives an instantaneous dose \(D=100\ \mathrm{mg}\) at time zero, with \(A(0^+)=F D\), then follows \(dA/dt=-kA\). Take \(F=0.80\), \(k=0.30\ \mathrm{h}^{-1}\), and calculate the amount and sensitivity to \(F\) at 2.00 h.
Differentiate the initial condition. The post-dose initial amount depends on \(F\), so \(S_F(0^+)=\partial(FD)/\partial F=D=100\ \mathrm{mg}\) per unit \(F\). After the dose, \(f=-kA\) has no direct dependence on \(F\), and the sensitivity satisfies \(dS_F/dt=-kS_F\). Therefore \(S_F(t)=D e^{-kt}\).
Calculate both quantities. Since \(e^{-0.30(2.00)}=e^{-0.6}=0.5488116\), the amount is \(A(2)=F D e^{-kt}=0.80(100)(0.5488116)=43.9049\ \mathrm{mg}\). The sensitivity is \(S_F(2)=100(0.5488116)=54.8812\ \mathrm{mg}\) per unit \(F\).
The dose affects the sensitivity through the parameter-dependent initial condition even though \(F\) is absent from the post-dose elimination equation. Setting the sensitivity initially to zero would therefore miss the dose's entire effect.
Worked Example: Exchange-Rate Sensitivity in Two Compartments
For a hypothetical two-compartment amount model, let \(A_1\) and \(A_2\) be amounts in mg and let \(dA_1/dt=-(k_{10}+k_{12})A_1+k_{21}A_2\) and \(dA_2/dt=k_{12}A_1-k_{21}A_2\). Use \(k_{10}=0.10\ \mathrm{h}^{-1}\), \(k_{12}=0.20\ \mathrm{h}^{-1}\), \(k_{21}=0.10\ \mathrm{h}^{-1}\), and initial amounts \(A_1(0)=100\ \mathrm{mg}\), \(A_2(0)=0\). Find a forward-Euler approximation after one step of \(\Delta t=0.10\) h for the states and their sensitivities to \(k_{12}\).
Set up the sensitivity system. With \(S_1=\partial A_1/\partial k_{12}\) and \(S_2=\partial A_2/\partial k_{12}\), the state Jacobian is \(\begin{pmatrix}-(k_{10}+k_{12})&k_{21}\\ k_{12}&-k_{21}\end{pmatrix}\), and the direct derivative is \(\begin{pmatrix}-A_1\\ A_1\end{pmatrix}\). Thus \(\dot S=J S+(-A_1,A_1)^\mathsf{T}\). Both initial sensitivities are zero because the stated initial amounts do not depend on \(k_{12}\).
Evaluate at time zero. The state rates are \(\dot A_1=-(0.30)(100)+(0.10)(0)=-30\ \mathrm{mg/h}\) and \(\dot A_2=(0.20)(100)-(0.10)(0)=20\ \mathrm{mg/h}\). Since \(S(0)=0\), the sensitivity rates are \(\dot S_1=-100\ \mathrm{mg}\) and \(\dot S_2=100\ \mathrm{mg}\), where the sensitivity-rate units are mg per unit \(k_{12}\) per hour.
Take one Euler step. \(A_1(0.10)\approx100+0.10(-30)=97\ \mathrm{mg}\), and \(A_2(0.10)\approx0+0.10(20)=2\ \mathrm{mg}\). Likewise, \(S_1(0.10)\approx0+0.10(-100)=-10\ \mathrm{mg\,h}\), and \(S_2(0.10)\approx0+0.10(100)=10\ \mathrm{mg\,h}\). These are first-order numerical approximations, not exact values. The opposite sensitivity signs are consistent with increasing \(k_{12}\) shifting amount out of compartment 1 and into compartment 2 over this initial short interval.
From State Sensitivities to PBPK Outputs
A model output may be a concentration rather than an amount. If \(C=A/V\), where both amount \(A\) and volume \(V\) may depend on \(\theta\), differentiating the ratio gives \(\partial C/\partial\theta=S_A/V-A(\partial V/\partial\theta)/V^2\). If the volume is fixed with respect to the parameter, the second term is zero. This distinction matters when a parameter changes both a state trajectory and a physiological quantity used to calculate the reported output.
For a general output \(y=h(t,x,\theta)\), its derivative is \(dy/d\theta=h_xS+h_\theta\). Thus the state sensitivity alone is not necessarily the sensitivity of a concentration, exposure summary, or other derived output. The output's explicit parameter dependence must also be differentiated.
Numerical Checks, Alternatives, and Limitations
Forward sensitivities are numerical solutions, so they require the same care as the state equations. The state and sensitivity scales can differ substantially; tolerances suitable for amounts may not resolve small derivatives. Check that tightening solver tolerances changes the sensitivities negligibly for the intended use. Verify the state trajectory separately, and inspect whether each dosing event was applied exactly once at the correct time.
A useful diagnostic is to compare selected forward derivatives with finite differences over a range of small perturbations. Agreement over a stable range supports the implementation; disagreement may indicate an incorrect derivative, omitted dose or initial-condition contribution, insufficient numerical accuracy, or a perturbation large enough to leave the local regime. Finite differences remain a practical alternative when analytic derivatives are difficult to implement, but can be sensitive to perturbation size and solver error. Forward equations avoid choosing a difference step, at the cost of deriving and integrating additional equations.
The method assumes differentiable dynamics between events and a consistent event structure in the parameter neighborhood of interest. Threshold-triggered dosing, discontinuous switches, or parameters that change which events occur can make ordinary local derivatives undefined or require event-specific treatment. Forward sensitivities also scale with the number of parameters: each parameter needs a sensitivity vector for all states. That computational cost can be substantial in a large PBPK system.
Practical Implementation Sequence
Specify the PBPK states, parameters, continuous inputs, initial conditions, and fixed dosing events.
For each target parameter, calculate the state Jacobian, direct parameter derivative, and initial sensitivity.
Integrate the augmented state-sensitivity system between events; update both states and sensitivities at each dose.
Check solver stability, compare selected results with finite differences, and differentiate any derived output formula.
Check Your Understanding
Answer each question, then compare with the solutions below.
- If the initial state is independent of \(\theta\), what is \(S(0)\)?
- For \(dA/dt=-kA+R(\theta)\), what direct input term appears in the sensitivity equation for \(R\)?
- A fixed-time bolus has \(A^+=A^-+FD\). What is the sensitivity jump with respect to \(F\)?
- Why can the sensitivity of a concentration differ from the sensitivity of its amount?
Solutions. (1) \(S(0)=0\), because the derivative of a parameter-independent initial state is zero. (2) The direct term is \(dR/d\theta\), so \(\dot S=-kS+dR/d\theta\) when \(k\) is fixed. (3) \(S_F^+=S_F^-+D\), since differentiating \(FD\) with respect to \(F\) gives \(D\). (4) If \(C=A/V\), then \(dC/d\theta=S_A/V-A(dV/d\theta)/V^2\); a parameter-dependent volume contributes an additional term.