Why Compare These Control Architectures?
A cytokine time course can rise and then decline under more than one biological explanation. The decline might follow reduced upstream stimulation, or it might arise because the response induces an inhibitor that later suppresses cytokine production. These are different causal structures, even when both can reproduce the observed concentrations. A useful comparison therefore asks not only which model fits the measured course, but also whether its residuals are credible and whether the models predict different outcomes under a targeted perturbation.
Here, \(X(t)\) denotes a measured activating cytokine or upstream signal, \(Y(t)\) the cytokine response of interest, and \(I(t)\) a candidate inhibitory mediator. Concentrations are in pg/mL and time is in hours. The examples use entirely hypothetical data and fitted outputs.
Two Explicit Model Hypotheses
A minimal feedforward model and a corresponding feedback model are:
The parameters \(\beta\), \(k_Y\), \(\alpha\), \(\delta_I\), and \(K_I\) are positive. The feedback term makes the production rate of \(Y\) decrease as \(I\) increases; \(K_I\) is the inhibitor concentration at which the production term is halved, holding \(X\) fixed. Both models require initial conditions. In the comparison below, the starting value of \(Y\) is measured, and the feedback model's starting value of \(I\) is fixed from baseline measurement or specified in advance. Treating an unknown initial value as another fitted parameter must be counted as additional flexibility.
For a fair calibration comparison, both models use the same \(X(t)\) observations, the same \(Y\) observations, the same time points, and the same observation-error assumption. The example assumes independent measurement errors with a common standard deviation of 0.5 pg/mL. This is plausible for the hypothetical data because each time point represents measurements from separate replicate wells. With repeated measurements from the same subject or culture, independence may not be plausible.
Fitting and Checking the Shared Time Course
Let \(y_j\) be the observed concentration of \(Y\) at time \(t_j\), and let \(\hat y_j\) be the model-predicted concentration at that time after numerically solving its ODEs. Define the residual as observed minus predicted concentration. Under the stated Gaussian error assumption, the residual contribution to the likelihood is proportional to the following weighted sum of squares:
This objective follows directly from the independent normal observation model. The density for one measurement is proportional to \(\exp[-(y_j-\hat y_j)^2/(2\sigma^2)]\). Multiplying the densities across observations adds the squared terms in the exponent, so maximizing the likelihood is equivalent to minimizing \(Q\) when \(\sigma\) is fixed. If the error scale is estimated rather than fixed, the fit still seeks small residuals, but the model comparison must account for the estimated scale and for differences in parameter flexibility.
Worked Example: Fitting Both Models to the Same Measurements
Suppose the measured \(Y\) concentrations at six times are as shown below. Numerical ODE integration and optimization give the displayed illustrative predictions for each model. They are fit outputs for this hypothetical data set, not algebraic solutions of the ODEs.
| Time (h) | Observed \(Y\) (pg/mL) | Feedforward prediction | Feedback prediction | Feedback residual (pg/mL) |
|---|---|---|---|---|
| 0 | 0.4 | 0.5 | 0.4 | 0.0 |
| 2 | 2.7 | 2.4 | 2.6 | 0.1 |
| 4 | 5.5 | 4.9 | 5.3 | 0.2 |
| 6 | 6.2 | 6.0 | 6.3 | -0.1 |
| 8 | 4.8 | 5.1 | 4.7 | 0.1 |
| 10 | 3.2 | 3.8 | 3.4 | -0.2 |
For the feedforward fit, the residuals are \(-0.1, 0.3, 0.6, 0.2, -0.3,\) and \(-0.6\) pg/mL. The unweighted residual sum of squares is
For the feedback fit, the residuals are \(0.0, 0.1, 0.2, -0.1, 0.1,\) and \(-0.2\) pg/mL, giving
With \(\sigma=0.5\) pg/mL, the weighted sums are \(Q_{\mathrm{FF}}=0.95/0.25=3.80\) and \(Q_{\mathrm{FB}}=0.11/0.25=0.440\). Thus, in this example, the feedback model has the closer calibration fit. That result alone does not establish that feedback is the true biological mechanism: the feedback model has more parameters and can be more flexible. The improvement should be considered alongside a complexity-aware model comparison and the models' predictive performance under conditions not used for calibration.
Residuals should also be inspected against time, fitted concentration, and relevant inputs. A good fit does not require residuals to equal zero, but they should not show systematic structure unexplained by the assumed error model. Here, the feedforward residuals are positive around 4 hours and negative at the last two times, suggesting that its fitted trajectory misses part of the time-course shape. The feedback residuals are smaller and alternate in sign, with no obvious trend in this small example. Six observations are too few to establish that the residual assumptions are correct.
What the Fit Can—and Cannot—Tell You
The calibration comparison is informative only under its assumptions: the specified ODE structure adequately represents the relevant biology; \(X(t)\) is measured well enough to use as an input; the initial conditions and units are correct; the numerical solver is accurate; and the observation-error model is reasonable. If assay variability grows with concentration, a common additive standard deviation may misweight the data. If measurements are censored, correlated, or collected from the same experimental unit over time, the residual model must reflect that design rather than treating all points as independent with constant variance.
A small residual sum of squares can conceal structural failure. Plot residuals over time: long runs of the same sign can indicate a missed delay, turnover process, or changing input. Plot residuals against predictions: increasing spread can indicate concentration-dependent error. Check that predictions remain biologically feasible and that fitted parameters are not stuck at arbitrary bounds. Also check numerical integration tolerances; a solver error comparable to the observation error can distort both the fit and the residual pattern.
Parameter identifiability is especially important in the feedback model. Over a short observation window, a strong inhibitor with rapid removal can resemble a weaker inhibitor with slower removal. If several parameter combinations produce nearly the same \(Y(t)\), the measured time course may identify the feedback-shaped trajectory without separately identifying \(\alpha\), \(\delta_I\), and \(K_I\). Measuring \(I\), adding informative time points, or perturbing the inhibitor can help distinguish these possibilities. A flexible model that fits better may still provide uncertain or unstable mechanistic parameters.
The main alternatives are not interchangeable. A calibration-only comparison measures agreement with the observed course, but does not directly test an intervention prediction. A complexity-aware comparison guards against rewarding extra adjustable parameters, but still depends on the candidate model set and error assumptions. A targeted perturbation tests a mechanistic consequence and can be more discriminating, although its interpretation depends on the perturbation's selectivity and on how well the intervention is represented in each model. These approaches complement one another.
Blockade Predictions as a Mechanistic Test
Consider blocking \(I\), the proposed inhibitory mediator, while holding the activating input \(X(t)\) fixed. In the feedback model, ideal complete blockade is represented by removing the inhibitory effect from the \(Y\)-production term, equivalent to setting \(I/K_I=0\) during blockade. In the feedforward model, which has no \(I\)-mediated path, the same blockade has no direct effect. This is a deliberately sharp contrast; in a real system, indirect effects or off-target activity could make the response differ.
Worked Example: Predicting the Response to Inhibitor Blockade
Suppose both fitted models predict a control concentration near 4.0 pg/mL at 8 hours. Under complete \(I\) blockade, the numerical predictions are 3.9 pg/mL for the feedforward model and 6.1 pg/mL for the feedback model. The predicted treatment effects, defined as blockade minus control, are therefore:
The model predictions differ by \(6.1-3.9=2.2\) pg/mL under blockade. If the standard deviation of a single measurement is 0.8 pg/mL, that separation is \(2.2/0.8=2.75\) measurement standard deviations. It is a substantial predicted contrast relative to the assumed measurement noise, not a guarantee that a study will detect the difference. Biological variability, imperfect blockade, uncertainty in fitted parameters, and uncertainty in the model structures can all reduce the effective separation.
The interpretation is conditional: if a selective blockade produces a clear increase in \(Y\) relative to control, that result is more consistent with the proposed negative-feedback path than with the specified feedforward model. It does not prove that \(I\) is the only regulator, nor does it establish causation unless the intervention and comparison support that causal interpretation.
Choosing a Discriminating Follow-Up
A useful follow-up should target a prediction on which the models disagree, not merely collect more observations where their trajectories overlap. For a simple planning calculation, let \(\mu_{\mathrm{FB},j}\) and \(\mu_{\mathrm{FF},j}\) be the two models' predicted mean concentrations at candidate time \(t_j\) under the same intervention. If measurement errors at the selected times are independent with known standard deviation \(\sigma\), a simple separation score is the sum of squared standardized differences:
This score expresses how far apart the predicted means are relative to measurement noise. It is a planning aid, not a complete power calculation or a substitute for propagating parameter and model uncertainty. Correlated measurements require a covariance-aware calculation rather than the independent-error formula above.
Worked Example: Selecting a Sampling Time Under Blockade
Suppose the predicted blockade concentrations are as follows, with a measurement standard deviation of 0.8 pg/mL at each time.
| Time (h) | Feedforward prediction (pg/mL) | Feedback prediction (pg/mL) | Difference (pg/mL) |
|---|---|---|---|
| 4 | 5.0 | 5.4 | 0.4 |
| 8 | 3.9 | 6.1 | 2.2 |
| 12 | 2.1 | 3.3 | 1.2 |
For these three times, the separation score is
The 8-hour sample contributes most of this score because its predicted difference is largest relative to the noise. If resources allow only one sampling time and the stated assumptions are credible, 8 hours is the strongest of these candidates. A more robust design might include an earlier sample to confirm comparable starting responses and an 8-hour sample for discrimination. Before committing, check whether blockade is expected to be established by 8 hours, whether the assay is reliable in that concentration range, and whether parameter uncertainty could make the two predictions overlap.
Common Pitfall and Practical Limits
Correction: A smaller residual sum of squares establishes only a closer fit to these observations under the chosen error model. It does not account by itself for the feedback model's extra flexibility, rule out omitted mechanisms, or demonstrate the intervention response. Compare models with appropriate complexity control, inspect residuals, and seek evidence from a perturbation that makes their predictions diverge.
Blockade experiments also have practical limits. A nominally selective antibody may not fully block its target, may alter other pathways, or may change the upstream input \(X(t)\). Measuring \(X\) and, where feasible, \(I\) alongside \(Y\) can reveal whether the intervention behaved as assumed. If the follow-up response falls between the two predictions, that is not automatically a failed experiment: it may indicate partial blockade, a hybrid network, or limitations in both candidate models. The next step is to diagnose which assumptions failed rather than forcing the result into one of the original alternatives.
Check Your Understanding
Use the stated model definitions and assumptions to answer each question. Brief solutions follow the questions.
- In the feedback equations, what is the role of \(I\), and what is the expected direct effect of complete \(I\) blockade on \(Y\) production?
- For residuals \(-0.2\) and \(0.3\) pg/mL with \(\sigma=0.5\) pg/mL, calculate their combined contribution to \(Q\).
- At a candidate time, the models predict 4.5 and 5.7 pg/mL, and the measurement standard deviation is 0.6 pg/mL. Calculate the standardized separation.
- Why is a lower calibration SSE insufficient to establish the feedback mechanism?
Solutions. (1) \(I\) is induced by \(Y\) and suppresses its production; complete blockade removes that modeled suppression, so the direct predicted effect is increased \(Y\) production relative to the unblocked feedback case. (2) The contribution is \((-0.2/0.5)^2+(0.3/0.5)^2=0.16+0.36=0.52\). (3) The predicted difference is \(5.7-4.5=1.2\) pg/mL, so the separation is \(1.2/0.6=2.0\) measurement standard deviations. (4) SSE alone does not penalize additional flexibility, verify residual assumptions, rule out alternative structures, or test a mechanism-specific perturbation prediction.