Diagnosing a Weak Pathway Rate Combination
A calibration can fit the observed pathway outputs closely and still leave some biological parameters poorly determined. The key question is not just whether the predictions match, but which changes in the parameters can occur with little change in the predictions. Fisher information describes local sensitivity to such changes; profile likelihood shows how the fit behaves when a parameter or parameter combination is fixed and the others are re-optimized.
The synthetic example below is deliberately transparent: its predictions and profiles can be calculated exactly. It illustrates a common pathway-model issue: data may determine the total flow through two branches much more precisely than how that flow is divided between the branches. The model is not a real study or a claim about a particular biological system.
Model, Data, and Assumptions
Let \(A\) be a normalized upstream pool held at \(A(t)=1\), feeding two downstream branches. The branch states \(B\) and \(C\) follow
The unknown rates \(k_1\) and \(k_2\) are nonnegative and have units of day\(^{-1}\). At \(t=1\) day, the assay reports the sum and a low-precision contrast:
The factor \(0.1\) represents a weakly scaled contrast measurement. The synthetic observations are \(S_{\mathrm{obs}}=1.000\) and \(D_{\mathrm{obs}}=-0.0200\), in normalized assay units. Assume independent Gaussian observation errors with known standard deviations \(\sigma_S=0.0200\) and \(\sigma_D=0.100\), respectively. The observation model, conditional on the rates, is therefore \(S_{\mathrm{obs}}\sim N(S,\sigma_S^2)\) and \(D_{\mathrm{obs}}\sim N(D,\sigma_D^2)\). There are no additional random effects or unknown error parameters in this example.
Worked Example: Solving the ODE and Finding the Fit
Because \(A(t)=1\), integrating each state equation from zero to \(t\) gives \(B(t)=k_1t\) and \(C(t)=k_2t\). At one day, the predictions are \(S=k_1+k_2\) and \(D=0.1(k_1-k_2)\).
The Gaussian negative log likelihood, after dropping terms constant in the rates, is one-half the weighted residual sum of squares:
Both squared terms are nonnegative. They are simultaneously zero when \(k_1+k_2=1.000\) and \(k_1-k_2=-0.200\). Solving gives \(\hat{k}_1=0.400\) day\(^{-1}\) and \(\hat{k}_2=0.600\) day\(^{-1}\), which satisfy the nonnegative-rate constraints. The predicted values are \(S=1.000\) and \(D=-0.0200\), exactly matching the hypothetical observations. An exact fit here is a consequence of two observations and two freely estimated rates, not evidence that the rates are precisely known.
Scaled Fisher Information
The unscaled sensitivity matrix uses derivatives with respect to the rates in day\(^{-1}\). Since the predictions are linear in the rates, these derivatives are constant throughout the parameter space. Scaling the columns by the fitted rate values makes each column correspond to a relative rate change: define \(\eta_j=(k_j-\hat{k}_j)/\hat{k}_j\) locally, and let \(D_k=\operatorname{diag}(\hat{k}_1,\hat{k}_2)\). The scaled information is \(I_{\mathrm{sc}}=D_k I D_k\). This is a local relative-change scaling, not a claim that the rates themselves are log-transformed.
Worked Example: Calculating and Interpreting the Scaled Information
The prediction derivatives with respect to \(k_1,k_2\) are \((1,1)\) for \(S\) and \((0.1,-0.1)\) for \(D\). Dividing each row by its observation standard deviation gives
Thus \(I=J^\mathsf{T}J\), and scaling by the fitted rates gives
The eigenvalues of the scaled matrix are approximately \(1300.08\) and \(0.4431\), a ratio of about 2934. The eigenvector for the smaller eigenvalue is approximately \((0.832,-0.555)\) in relative-change coordinates. Multiplying its components by the fitted rates gives approximately \((0.333,-0.333)\) day\(^{-1}\): an increase in one rate offset by a decrease in the other. This direction changes the branch allocation while leaving the total rate nearly unchanged. Its small information eigenvalue signals weak local curvature, not that either individual rate has no biological meaning.
For this known-variance Gaussian model, \(I=J^\mathsf{T}J\) is the expected Fisher information and also the curvature matrix of the weighted least-squares objective. The eigenvalues depend on parameter scaling, so reporting them without stating the scaling is incomplete. The eigenvectors also need interpretation in the model's parameter coordinates: here the weak direction corresponds to an opposing change in the two rates, not an equal change.
Profiles Show the Shape of the Weak Direction
A profile likelihood fixes a quantity of interest and maximizes the likelihood over the remaining nuisance parameters. In this example, write \(Q=-2\{\ell-\ell(\hat{k}_1,\hat{k}_2)\}\), the increase in negative twice the log likelihood from its optimum. A profile for \(k_1\) is obtained by fixing \(k_1\) and minimizing \(Q\) over \(k_2\geq0\). A profile for a derived quantity such as the rate difference is constructed in the same way, fixing that quantity while allowing the other free direction to adjust.
Worked Example: Profiling the First Rate
Let \(a=k_1-0.400\) and \(b=k_2-0.600\). Then the objective difference is
At fixed \(a\), differentiation with respect to \(b\) gives \(5000(a+b)+2(a-b)=0\), so the unconstrained nuisance optimum is \(b^*=-2499a/2501\). Substitution gives
For example, at \(k_1=0.800\) day\(^{-1}\), \(a=0.400\) day\(^{-1}\). The optimizing nuisance value is \(k_2=0.600-(2499/2501)(0.400)\approx0.2003\) day\(^{-1}\), which is feasible. The profile increase is \(Q_{\mathrm{prof}}\approx3.9984(0.400)^2=0.6397\). The increase is small despite the first rate moving by half its fitted value; the second rate compensates to preserve the well-measured sum.
The following values are calculated from this profile formula. The nuisance optimum is nonnegative at each listed point, so the nonnegative-rate constraint does not alter these entries.
| Fixed \(k_1\) (day\(^{-1}\)) | Optimizing \(k_2\) (day\(^{-1}\)) | Profile increase \(Q\) |
|---|---|---|
| 0.200 | 0.7998 | 0.1599 |
| 0.400 | 0.6000 | 0.0000 |
| 0.600 | 0.4002 | 0.1599 |
| 0.800 | 0.2003 | 0.6397 |
| 1.000 | 0.000480 | 1.439 |
Plotting the fixed values on the horizontal axis and \(Q_{\mathrm{prof}}\) on the vertical axis gives a shallow, approximately quadratic curve near the optimum. For comparison, define the sum \(u=k_1+k_2\) and difference \(v=k_1-k_2\). The objective separates exactly:
When \(v\) is fixed, the optimum is \(u=1.000\), provided the implied rates are nonnegative. Over the feasible range \(-1\leq v\leq1\) at this sum, the difference profile is \(Q_{\mathrm{prof},v}=(v+0.200)^2\). Its values at \(v=-1.0,-0.6,-0.2,0.2,0.6,1.0\) are \(0.640,0.160,0,0.160,0.640,1.440\). The profile for \(u\), in contrast, is \(2500(u-1.000)^2\), which rises sharply: a change of \(0.0200\) in the sum already gives \(Q=1\). Together, these profiles make the weak contrast and strong sum visually and quantitatively clear.
How to Diagnose the Calculation
A reliable diagnosis checks the model evaluation, sensitivities, scaling, and profiles as related but distinct calculations. A broad profile can be genuine; an optimizer that fails to re-optimize nuisance parameters can also create a misleading profile. In a numerical ODE model, practical checks include:
- Verify predictions at the fitted rates. Re-evaluate the ODE and observation model, and confirm that the residuals and objective agree with the reported fit.
- Check sensitivities independently. Compare analytic or sensitivity-equation derivatives with finite differences over more than one perturbation size. Use perturbations appropriate to the parameter scale, and check that solver tolerances are sufficiently tight.
- Check matrix construction. Standardize each sensitivity row by the observation standard deviation before forming \(J^\mathsf{T}J\). Confirm that a stated parameter scaling is applied on both sides of the information matrix.
- Check each profile optimization. Start nuisance optimization from multiple plausible values, inspect convergence and constraints, and verify the objective at the resulting nuisance estimates. A profile that changes materially with starting values needs numerical investigation before biological interpretation.
- Vary the profile grid and range. A coarse grid can conceal a narrow feature, while an overly narrow range can make a flat direction look well determined simply because it was not explored.
For the exact linear example, derivatives and profiles can be checked algebraically. In nonlinear pathway models, local Fisher information summarizes curvature near one fitted point; profiles can reveal asymmetry, changing curvature, or constraint effects farther away. This is a comparison of complementary diagnostics, not a reason to expect every profile to be quadratic.
Interpretation, Alternatives, and Limits
The estimate of the sum is informed chiefly by the precise sum measurement, while the estimate of the contrast depends on the much less precise contrast assay. The weak direction is therefore a feature of the observation design and error assumptions, not necessarily a defect in the ODE solver. An experiment that measured the two branches separately with adequate precision could add information about their allocation; simply collecting more precise measurements of the same sum would mainly reinforce the already strong direction.
The local Fisher calculation is compact and useful for identifying sensitive combinations near a fit. Profile likelihood is more computationally demanding, but tests how the likelihood changes as a chosen parameter or combination moves while other parameters adjust. In this constructed model, the quadratic objective makes the profile results exact. In a nonlinear model, profiles can be asymmetric or affected by parameter boundaries, and the local eigenvectors need not describe the whole profile. Neither diagnostic alone proves that a parameter is biologically identifiable under all possible experiments.
A likelihood-ratio cutoff such as \(Q=3.84\) is often used for an approximate 95% one-parameter interval under regular large-sample conditions. Here there are only two synthetic observations, and nonnegative parameter constraints can affect profile behavior. Therefore that cutoff should not be treated as automatically calibrated for this small example. The calculations diagnose relative curvature and compensation; they do not establish a calibrated confidence interval.
The design of this example also imposes an important limitation: there are exactly two independent observations and two rates, with known error standard deviations. The conclusions are conditional on that observation model. Additional unknown error parameters, correlated errors, uncertain upstream input, or a different measurement scale would change the likelihood and its information. The formulas must be rebuilt from the actual model and observation assumptions rather than transferred unchanged.
Check Your Understanding
Answer each question using the model and assumptions above. Short solutions are included after the questions.
- What rate combination is most precisely estimated, and which is weakly estimated?
- At fixed \(k_1=0.600\) day\(^{-1}\), find the profiled value of \(k_2\) and the increase \(Q\).
- Why must the sensitivities be divided by their observation standard deviations before constructing the stated Fisher information?
- If the measurement of the sum became more precise but the contrast assay stayed unchanged, which combination would benefit most directly?
- Why can a shallow profile not be attributed to biological weak estimation until nuisance optimization has been checked?
Solutions. (1) The sum \(u=k_1+k_2\) is most precisely estimated; the difference \(v=k_1-k_2\), corresponding to opposing rate changes, is weak. (2) Here \(a=0.200\), so \(k_2=0.600-(2499/2501)(0.200)\approx0.4002\) day\(^{-1}\), and \(Q=(10000/2501)(0.200)^2\approx0.1599\). (3) Standardization weights residual changes by their Gaussian variances, so an equal raw prediction change carries less information when its observation error is larger. (4) The sum would benefit most directly because the precise sum observation contributes information along \(u\). (5) Failed or incomplete nuisance re-optimization can leave the objective artificially high or low at fixed parameter values; convergence, multiple starts, and direct objective checks are needed to distinguish that numerical problem from a genuinely flat likelihood direction.