Why Finite Differences Need Care
A local sensitivity is a derivative of a model output with respect to a parameter at a specified parameter set. In practice, a PBPK model may not provide that derivative directly. One simple approach is to perturb one parameter, rerun the model, and compare the resulting outputs. The calculation is straightforward; choosing a perturbation that produces a trustworthy estimate is not.
Finite differences are useful when the model can be evaluated but a derivative is unavailable. They can also serve as a check on other derivative calculations. However, a very large perturbation averages behavior over a wider parameter range, while a very small perturbation can make solver error dominate the difference. The goal is therefore not to choose the smallest possible step, but to find a feasible step for which the estimate is both local and numerically stable.
This tutorial concerns parameter perturbations, not the discretization of time in the PBPK differential equations. The output \(g\) might be a concentration at a specified time, an exposure summary, or another defined endpoint. The endpoint, dosing scenario, and values of all unperturbed parameters must be held fixed across evaluations. As in the earlier material on dimensionless sensitivity coefficients, the derivative is local; unlike a dimensionless coefficient, it retains the units of output divided by parameter.
Three Difference Formulas and Their Assumptions
For a scalar parameter \(\theta\), write the model output as \(f(\theta)\), with all other inputs fixed. Let \(h>0\) be a perturbation in the same units as \(\theta\). The forward difference uses the baseline and an increased parameter value; the backward difference uses the baseline and a decreased value. The central difference uses both directions.
These estimates require the relevant perturbed values to be within the model’s valid parameter domain. Their usual error-order statements also require smoothness near \(\theta\). If \(f\) has a locally bounded third derivative, Taylor expansion gives the forward and backward estimates as the true derivative plus a first-order error in \(h\), with a remainder of order \(h^2\). The central estimate has a second-order error. If only \(f\in C^2\) is known, the forward and backward remainders after the second-derivative terms are generally \(o(h)\), not necessarily \(O(h^2)\); do not claim an \(O(h^2)\) remainder without stronger regularity.
If \(f\) has a locally bounded third derivative, Taylor’s theorem gives \(f(\theta+h)=f(\theta)+hf'(\theta)+\tfrac{h^2}{2}f''(\theta)+O(h^3)\) and \(f(\theta-h)=f(\theta)-hf'(\theta)+\tfrac{h^2}{2}f''(\theta)+O(h^3)\).
Subtracting the baseline and dividing by \(h\) yields \(D_{\mathrm{forward}}(h)=f'(\theta)+\tfrac{h}{2}f''(\theta)+O(h^2)\) and \(D_{\mathrm{backward}}(h)=f'(\theta)-\tfrac{h}{2}f''(\theta)+O(h^2)\). The leading errors have opposite signs.
Subtract the two Taylor expansions and divide by \(2h\). The terms involving \(f''(\theta)\) cancel, leaving \(D_{\mathrm{central}}(h)=f'(\theta)+O(h^2)\). This cancellation explains central differences’ higher accuracy for a smooth output, not immunity to numerical error.
More precisely, if \(|f'''(u)|\leq M_3\) throughout the interval from \(\theta-h\) to \(\theta+h\), the magnitude of the central truncation error is bounded by \(M_3h^2/6\). For the forward or backward difference, a bound \(|f''(u)|\leq M_2\) gives a first-order truncation-error bound \(M_2h/2\). These are truncation errors for exact evaluations of \(f\), not total errors in a computational PBPK result.
Worked Example: Comparing the Three Formulas
Suppose a hypothetical PBPK endpoint is described near the parameter of interest by \(f(\theta)=12-0.8\theta+0.06\theta^2\), where \(\theta\) is clearance in L/h and \(f\) is AUC in mg·h/L. At \(\theta=5\) L/h, estimate the local derivative using \(h=0.5\) L/h. The function is smooth and all three parameter values are in its domain.
The baseline output is \(f(5)=12-0.8(5)+0.06(25)=9.5\) mg·h/L. The perturbed outputs are \(f(5.5)=12-0.8(5.5)+0.06(5.5^2)=9.415\) mg·h/L and \(f(4.5)=12-0.8(4.5)+0.06(4.5^2)=9.615\) mg·h/L.
The forward estimate is \((9.415-9.5)/0.5=-0.170\) (mg·h/L)/(L/h). The backward estimate is \((9.5-9.615)/0.5=-0.230\) (mg·h/L)/(L/h). The central estimate is \((9.415-9.615)/(2(0.5))=-0.200\) (mg·h/L)/(L/h).
Differentiating the specified function gives \(f'(\theta)=-0.8+0.12\theta\), so \(f'(5)=-0.200\) (mg·h/L)/(L/h). The central estimate is exact here because the output function is quadratic; the two one-sided estimates differ because each has a curvature-related first-order error. A general PBPK output need not be quadratic, so exactness in this example should not be assumed elsewhere.
Choosing a Perturbation Size
A computed PBPK output is not exact. Suppose each numerical evaluation differs from the ideal model output by at most \(\varepsilon\) output units. In a forward or backward difference, the difference of two such errors, divided by \(h\), can contribute as much as \(2\varepsilon/h\). In a central difference, the two perturbed evaluations can contribute at most \(2\varepsilon/(2h)=\varepsilon/h\). Thus shrinking \(h\) reduces truncation error but increases the potential effect of solver error. If solver error is not bounded, these expressions are guides to the trade-off, not guaranteed error bounds.
For smooth functions with local bounds \(M_2\) and \(M_3\), simple worst-case bounds for total error are \(M_2h/2+2\varepsilon/h\) for a forward difference and \(M_3h^2/6+\varepsilon/h\) for a central difference. Minimizing these approximate bounds gives \(h\approx\sqrt{4\varepsilon/M_2}\) for forward differences and \(h\approx(3\varepsilon/M_3)^{1/3}\) for central differences. These expressions clarify how error sources interact, but \(M_2\), \(M_3\), and the effective output error are seldom known accurately in an applied model. A small pilot sweep across several feasible step sizes is usually more informative than treating these approximations as a universal step-size rule.
The size of \(h\) must also be meaningful for the parameter. For example, a step of \(0.01\) has different implications for a parameter near \(0.02\) than for one near \(100\). Use the parameter’s units and scale to choose a relative perturbation that changes the input enough to exceed numerical noise, while remaining local to the point of interest. Report the actual perturbation and the units of the resulting derivative. When parameters have bounds, check feasibility before running the perturbed model.
Worked Example: Truncation Error Versus Solver Error
Consider the hypothetical endpoint \(f(\theta)=4+2\theta+0.5\theta^2+0.1\theta^3\), where \(\theta\) is a physiological parameter in specified units and \(f\) is AUC in mg·h/L. At \(\theta=2\), the exact derivative is \(f'(2)=2+2+0.3(2^2)=5.2\) (mg·h/L) per parameter unit. Suppose numerical output error is no larger than \(0.0005\) mg·h/L at each evaluation. The parameter values in the perturbation sweep are assumed feasible.
For a central difference with \(h=0.2\), \(f(2.2)=4+4.4+0.5(2.2^2)+0.1(2.2^3)=11.8848\) and \(f(1.8)=4+3.6+0.5(1.8^2)+0.1(1.8^3)=9.8032\). Therefore \(D_{\mathrm{central}}(0.2)=(11.8848-9.8032)/0.4=5.204\). The estimate differs from the exact derivative by \(0.004\), its central-difference truncation error.
For this function, \(f'''(\theta)=0.6\), so the central truncation bound is \(0.6(0.2^2)/6=0.004\), matching that error. The solver-error contribution is bounded by \(0.0005/0.2=0.0025\) derivative units. At \(h=0.02\), the truncation error falls to \(0.6(0.02^2)/6=0.00004\), but the solver-error contribution can rise to \(0.0005/0.02=0.025\). A much smaller step is not necessarily more accurate. The two-term bound suggests a step near \(h=(3(0.0005)/0.6)^{1/3}\approx0.136\) parameter units, but a solver-tolerance check and a step-size sweep are still needed in an actual computation.
Parameter Bounds, Diagnostics, and Alternatives
At an interior parameter value, a central difference is often attractive when both perturbations are feasible and model evaluations are sufficiently smooth. At a lower or upper bound, one of those perturbations may be invalid. A forward difference at a lower bound or a backward difference at an upper bound can estimate a one-sided derivative, provided the model is defined and sufficiently smooth on the feasible side. That one-sided derivative is not automatically the same as a two-sided derivative if the model has a boundary kink or is defined only on one side.
For a parameter close to a bound, a central step may be feasible in only one direction or may require an impractically small \(h\). Do not silently clip a perturbation to the bound and still divide by \(2h\): the actual changes would no longer match the central-difference formula. Instead, use an appropriate one-sided formula with its actual step, choose a smaller symmetric step if it remains numerically useful, or clearly define a reparameterization and interpret the derivative in that coordinate. Reparameterization changes the derivative’s scale through the chain rule.
A practical diagnostic is to calculate the estimate at several steps, such as \(h\), \(h/2\), and \(2h\), while keeping the baseline and model settings fixed. For a smooth output with sufficiently accurate evaluations, central estimates should often approach a stable value as \(h\) is reduced, before solver error takes over. One-sided estimates generally have larger truncation error at the same step. If estimates change erratically, check whether outputs are reproducible, whether the solver tolerances are adequate, and whether the parameter perturbations cross a threshold, trigger a different event, or cause clipping or a change in model logic. Also confirm that all non-perturbed parameters, dosing inputs, and endpoint definitions are identical.
Finite differences are a general, black-box option: they need only model evaluations, but can require multiple runs per parameter and can be sensitive to numerical noise. Analytic derivatives, when available and correctly implemented, avoid differencing error but require derivative expressions. Automatic differentiation can propagate derivatives through supported calculations, but may not handle every solver or model component. Forward sensitivity equations and adjoint sensitivities are other derivative-calculation approaches; their detailed formulations are treated separately. Regardless of approach, a finite-difference comparison can be a useful diagnostic, but agreement is meaningful only when the perturbation and solver errors have been considered.
Worked Example at a Parameter Boundary
Worked Example: A One-Sided Derivative at a Lower Bound
Suppose a hypothetical AUC endpoint is \(f(\theta)=8-2\theta+\theta^2\) mg·h/L, where \(\theta\) is a dimensionless fraction constrained to \(0\leq\theta\leq1\). Estimate the derivative at the lower bound \(\theta=0\). The model is smooth on the feasible interval, and \(h=0.1\) remains within that interval.
A backward or central difference would require a value below zero, which is outside the stated parameter domain. The feasible forward estimate is \(D_{\mathrm{forward}}(0.1)=[f(0.1)-f(0)]/0.1\). Since \(f(0)=8\) and \(f(0.1)=8-2(0.1)+(0.1)^2=7.81\), the estimate is \((7.81-8)/0.1=-1.9\) (mg·h/L) per unit fraction.
The right-hand derivative is \(f'(0)=-2\) (mg·h/L) per unit fraction. The estimate’s difference of \(0.1\) is consistent with the forward formula’s curvature term: \(f''(\theta)=2\), so \(hf''(0)/2=0.1\). The estimate is a right-sided approximation at the boundary, not evidence that a negative parameter value is valid.
Check Your Understanding
Use the definitions and error reasoning above to answer each question. Complete solutions follow each problem.
- A smooth model output \(f(\theta)\) is evaluated with a central difference at an interior point. What is the leading truncation-error order if the third derivative is locally bounded, and what regularity warning applies if only \(f\in C^2\) is known?
Solution: With a locally bounded third derivative, the central difference has truncation error \(O(h^2)\). If only \(f\in C^2\) is known, the symmetric Taylor expansions imply an error \(o(h)\); that assumption alone does not justify claiming \(O(h^2)\).
- For \(f(\theta)=3+\theta^2\), calculate the forward difference at \(\theta=2\) using \(h=0.1\), and compare it with the derivative.
Solution: \(f(2)=7\) and \(f(2.1)=3+2.1^2=7.41\). Thus \(D_{\mathrm{forward}}(0.1)=(7.41-7)/0.1=4.1\). The exact derivative is \(f'(2)=2(2)=4\), so the estimate is \(0.1\) higher.
- If each model evaluation has absolute numerical error at most \(\varepsilon\), what is the maximum contribution of those errors to a central-difference estimate with step \(h\)?
Solution: Each of the two perturbed evaluations can err by up to \(\varepsilon\), so their difference can err by up to \(2\varepsilon\). Dividing by \(2h\) gives a maximum contribution of \(\varepsilon/h\).
- A parameter is at its lower bound and negative perturbations are invalid. Which basic difference is feasible, and what derivative does it approximate?
Solution: A forward difference is feasible if the positive perturbation remains within the domain. It approximates the right-hand derivative at the bound, with first-order truncation error under the usual smoothness conditions.