Sensitivity Matrices as Local Estimability Diagnostics

Why Examine a Sensitivity Matrix?

A fitted QSP model can reproduce observed data while leaving some parameters poorly determined. One reason is that, near the fitted parameter values, changing one parameter may barely change the predicted observations. Another is that changing two parameters in compensating directions may produce nearly the same output changes. The sensitivity matrix makes these local effects explicit: each column describes how the model’s observation predictions change when one parameter changes.

This is a local diagnostic, evaluated at a specified parameter vector and observation design. It does not by itself establish global structural identifiability or guarantee that a parameter will be estimated precisely from noisy data. It is most informative when its rows represent the actual observed outputs and its columns are scaled to meaningful parameter changes. The earlier material on Fisher information for ODE observation models provides the related statistical framework; here the emphasis is on interpreting the derivative matrix itself.

Definition: Let $\hat{\theta}=(\hat{\theta}_1,\ldots,\hat{\theta}_p)$ be a fitted parameter vector and let $g_i(\theta)$ be the model prediction for observation $i$, for $i=1,\ldots,m$. The observation-by-parameter sensitivity matrix at the fit is $$ J_{ij}=\left.\frac{\partial g_i(\theta)}{\partial \theta_j}\right|_{\theta=\hat{\theta}}, \qquad i=1,\ldots,m,\quad j=1,\ldots,p. $$ It has $m$ rows and $p$ columns. Entry $J_{ij}$ is the local change in predicted observation $i$ per unit change in parameter $\theta_j$, with all other parameters held fixed.

The units of a raw entry are the units of observation $i$ divided by the units of parameter $j$. Thus, raw column magnitudes cannot be compared fairly when parameters use different units or scales. A change of one hour$^{-1}$ may be enormous for one rate constant but negligible for another; likewise, a change of one concentration unit may be large or small depending on the output.

Scaling the Local Changes

Choose a positive scale $q_j$ for each parameter and a positive scale $r_i$ for each observation. A parameter scale might be a scientifically meaningful change, a typical parameter magnitude, or a standard-error scale. A response scale might be a characteristic output magnitude or, when assessing effects relative to measurement error, the observation’s standard deviation.

Definition: The scaled sensitivity matrix $J_s$ has entries $$ (J_s)_{ij}=\frac{q_j}{r_i}J_{ij}. $$ Each entry describes the change in the prediction, measured in units of $r_i$, caused locally by a change of $q_j$ in parameter $\theta_j$. Equivalently, for scaled parameter changes $u_j=(\theta_j-\hat{\theta}_j)/q_j$ and scaled prediction changes $v_i=(g_i(\theta)-g_i(\hat{\theta}))/r_i$, the first-order approximation is $v\approx J_su$.

The approximation follows from a first-order Taylor expansion at $\hat{\theta}$: $g_i(\hat{\theta}+\delta\theta)\approx g_i(\hat{\theta})+\sum_j J_{ij}\delta\theta_j$. Substituting $\delta\theta_j=q_j u_j$ and dividing the prediction change by $r_i$ gives $v\approx J_su$. This also clarifies what the diagnostic says: a column is the pattern of output changes produced by a one-scale-unit change in its parameter.

Scaling is a substantive choice, not a cosmetic adjustment. If parameters are positive and a relative change is appropriate, one can choose $q_j=\hat{\theta}_j$; this corresponds locally to a unit change in the parameter’s relative scale. That choice is not suitable when $\hat{\theta}_j=0$, and it may be inappropriate when estimates are near a boundary. In those cases, use a meaningful nonzero scale. If errors are independent with standard deviations $\sigma_i$, setting $r_i=\sigma_i$ expresses prediction changes relative to observation noise. For correlated errors, the observations can instead be scaled using the error covariance structure; a simple row-by-row division does not account for correlation.

Obtaining Sensitivities From an ODE Model

Suppose the state vector $x(t,\theta)$ solves $\dot{x}=f(x,\theta,t)$, with initial condition $x(0,\theta)=x_0(\theta)$. For a differentiable model and fixed observation times, differentiate the state equation with respect to parameter $\theta_j$. The state sensitivity $S_j(t)=\partial x(t,\theta)/\partial\theta_j$ satisfies

\dot{S}_j(t)=f_x(x,\theta,t)S_j(t)+f_{\theta_j}(x,\theta,t), \qquad S_j(0)=\frac{\partial x_0(\theta)}{\partial\theta_j}.

Here $f_x$ is the derivative of the vector field with respect to the state, and $f_{\theta_j}$ is its derivative with respect to parameter $\theta_j$; both are evaluated along the model solution. If observation $i$ is $g_i(\theta)=h(x(t_i,\theta),\theta,t_i)$, then

\frac{\partial g_i}{\partial\theta_j} = h_x(x(t_i),\theta,t_i)S_j(t_i) + h_{\theta_j}(x(t_i),\theta,t_i).

The second term accounts for a parameter entering the observation function directly, rather than only through the state. These equations assume the relevant functions are differentiable and that the observation times are fixed. Parameter-dependent event times, discontinuities, or dosing switches require the corresponding event or jump sensitivities; simply integrating the smooth equations across such events can give incorrect derivatives.

Worked Example: Scaling Sensitivities for Exponential Decay

Consider the hypothetical observation model $C(t)=A\exp(-kt)$, with fitted values $\hat{A}=10$ mg/L and $\hat{k}=0.5$ h$^{-1}$. Predictions are recorded at $t=0,1,2$ h. Use parameter scales $q_A=10$ mg/L and $q_k=0.5$ h$^{-1}$, and a response scale $r_i=10$ mg/L for each observation. The model and observation function are smooth, the times are fixed, and the scales are positive, so the local scaled derivative is defined.

The required derivatives are $\partial C/\partial A=\exp(-kt)$ and $\partial C/\partial k=-At\exp(-kt)$. Thus the scaled amplitude entry is $(10/10)\exp(-0.5t)$, and the scaled rate-constant entry is $(0.5/10)[-(10)t\exp(-0.5t)]=-0.5t\exp(-0.5t)$. Evaluating at the three times gives

J_s= \begin{pmatrix} 1 & 0\\ 0.6065 & -0.3033\\ 0.3679 & -0.3679 \end{pmatrix}.

For example, at $t=1$ h the rate-constant entry is $(0.5\ \mathrm{h}^{-1}/10\ \mathrm{mg/L})[-(10\ \mathrm{mg/L})(1\ \mathrm{h})\exp(-0.5)]=-0.3033$. At $t=2$ h it is $(0.5/10)[-(10)(2)\exp(-1)]=-0.3679$, with the units cancelling under the specified scales.

The scaled amplitude-column norm is $\sqrt{1^2+0.6065^2+0.3679^2}\approx1.226$; the scaled rate-column norm is $\sqrt{0^2+0.3033^2+0.3679^2}\approx0.4768$. Under these scales, a one-scale-unit change in $A$ has a larger total local effect across these predictions than a one-scale-unit change in $k$. The columns are not proportional, so their predicted patterns differ in this example. These calculations describe local response at the fitted values; they do not alone quantify uncertainty or prove that the parameters can be estimated reliably from noisy data.

Reading Columns and Checking the Diagnostic

A column with small entries under scientifically relevant scales indicates that the observations change little for a parameter-scale perturbation. That parameter’s effect may be weakly represented by this particular output and sampling design. If two columns are exactly proportional, their local prediction patterns are indistinguishable up to a rescaling: some small changes in one parameter can be offset by changes in the other. Near proportionality is a warning of potential local compensation, not proof that the fitted parameters are globally nonidentifiable.

The rank of the derivative matrix gives a basic algebraic check: if its columns are linearly dependent, at least one parameter-change direction produces no first-order change in the predicted observations. Full column rank rules out that particular exact first-order dependence, but it does not ensure useful precision. Small noise, poor scaling, nonlinear behavior away from the fit, or numerical error may still make inference unreliable. Column norms and comparisons between columns are therefore descriptive diagnostics, not stand-alone pass/fail criteria.

Worked Example: A Three-Parameter Response Model

Consider the hypothetical response $R(t)=B+E(1-\exp(-kt))$, fitted at $\hat{B}=5$, $\hat{E}=4$, and $\hat{k}=1$ day$^{-1}$. Take observations at $t=0,1,2$ days, parameter scales $q_B=5$, $q_E=4$, $q_k=1$ day$^{-1}$, and response scale $r_i=5$ response units. This smooth model has fixed observation times and nonzero scales, satisfying the conditions for the local derivative calculation.

The derivatives are $\partial R/\partial B=1$, $\partial R/\partial E=1-\exp(-kt)$, and $\partial R/\partial k=Et\exp(-kt)$. Applying the scale factors gives the rows at $t=0,1,2$:

J_s= \begin{pmatrix} 1 & 0 & 0\\ 1 & 0.5057 & 0.2943\\ 1 & 0.6917 & 0.2165 \end{pmatrix}.

For instance, at $t=1$, the scaled $k$ entry is $(1\ \mathrm{day}^{-1}/5)(4\ \mathrm{units})(1\ \mathrm{day})\exp(-1)=0.2943$. The displayed matrix has full column rank: subtracting its first row from each later row leaves a two-by-two determinant of $0.5057(0.2165)-0.6917(0.2943)\approx-0.0941$, which is nonzero. Thus these three predictions have distinct first-order effects for the three parameters in exact arithmetic. This rank calculation does not establish that those effects are large relative to measurement error or robust to numerical approximation.

Worked Example: Weak Rate Sensitivity at Early Times

Use the exponential-decay model and scales from the first example, but suppose observations are taken only at $t=0,0.1,0.2$ h. The model remains smooth, and the same positive scales apply. The scaled columns are calculated as $\exp(-0.5t)$ and $-0.5t\exp(-0.5t)$, producing

J_s\approx \begin{pmatrix} 1 & 0\\ 0.9512 & -0.04756\\ 0.9048 & -0.09048 \end{pmatrix}.

The amplitude-column norm is $\sqrt{1+0.9512^2+0.9048^2}\approx1.650$. The rate-column norm is $\sqrt{0.04756^2+0.09048^2}\approx0.1022$, about $0.0620$ times the amplitude-column norm. With these response and parameter scales, the early observations change much less for a scale-sized change in $k$ than for one in $A$. The zero rate sensitivity at time zero is expected: $C(0)=A$ does not depend on $k$. This is a local, design-dependent warning about the rate effect, not a universal claim that $k$ is unestimable under other sampling times or outputs.

Diagnostics, Alternatives, and Limitations

A reliable sensitivity calculation requires more than a matrix with plausible-looking entries. Check that the derivatives correspond to the fitted parameter vector, the correct observation function, and the actual sampling times. Verify units before comparing columns, and document the choices of $q_j$ and $r_i$. If sensitivity equations are integrated numerically, tighten solver tolerances and confirm that the entries stabilize. As an independent check, compare selected entries with finite differences: perturb $\theta_j$ by a small amount and compare the prediction change divided by that perturbation with $J_{ij}$. Finite differences can themselves be inaccurate if the step is too small (round-off) or too large (nonlinearity), so they are a check rather than an automatic gold standard.

The main alternatives answer related but different questions. Raw derivatives retain physical units and can be useful for a particular physical perturbation, but comparing raw column sizes across differently scaled parameters is misleading. Scaled derivatives provide a more interpretable local comparison, although their conclusions depend on the chosen scales. A noise-scaled or covariance-whitened derivative matrix emphasizes effects relative to observational uncertainty; it does not remove model misspecification or replace a full statistical analysis. Fisher information incorporates the observation-error model into a local information calculation, whereas the sensitivity matrix displays the underlying prediction derivatives directly.

This diagnostic is local: it describes a first-order approximation near one fitted point, so it may miss nonlinear behavior farther away. It can also depend strongly on the experiment—sampling times, doses, measured outputs, and perturbations determine which parameter effects are visible. Near a parameter boundary, a symmetric derivative may not describe feasible changes on both sides. Finally, a small column can reflect an unsuitable response scale as well as genuinely weak information. Sensitivity matrices guide scrutiny of a fit and its design; they do not replace checks of model assumptions, uncertainty, or predictions of scientific interest.

Common pitfall: A student sees a raw derivative column with small numbers and concludes that its parameter is poorly estimable. That conclusion ignores the parameter and response units: rescaling either can change the numerical column without changing the model. Choose and report meaningful $q_j$ and $r_i$, inspect scaled effects relative to observation noise when appropriate, and treat the result as a local diagnostic rather than a complete identifiability or uncertainty analysis.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions.

If a model has 12 predicted observations and 4 fitted parameters, what are the dimensions of its observation-by-parameter sensitivity matrix?
A raw derivative has units of concentration per day$^{-1}$. What units does its scaled entry have when multiplied by $q_j$ in day$^{-1}$ and divided by $r_i$ in concentration?
What does exact proportionality between two sensitivity columns imply locally, and what does it not prove?
Why might a zero sensitivity to a rate parameter at the first observation time be scientifically expected?
Give two checks that can reveal a numerical sensitivity calculation is unreliable.

Solutions. (1) It is $12\times4$, with one row per prediction and one column per parameter. (2) The units cancel: $(\text{concentration}/\text{day}^{-1})(\text{day}^{-1})/(\text{concentration})=1$, so the entry is dimensionless. (3) The two parameters have proportional first-order prediction patterns, permitting locally compensating changes; this alone does not prove global nonidentifiability or characterize behavior far from the fit. (4) For $C(t)=A\exp(-kt)$, $C(0)=A$, so $\partial C(0)/\partial k=0$. (5) For example, check whether entries stabilize under tighter ODE solver tolerances and compare selected entries against appropriately chosen finite-difference calculations.