Estimating Cascade Parameters From Response Time Courses

Why cascade fitting depends on the time course

A cascade model describes how a signal passes through successive stages before appearing in an observed response. Its parameters are not all learned equally from every observation. Early measurements help locate the rise and delay; measurements during the rising phase constrain the response shape; later measurements show whether the signal decays as predicted. If a time course covers only one of these regions, several combinations of rate constant, stage number, and gain can produce similar predictions.

Here the fitting target is a response following a known, sustained step input. For identical first-order stages, the normalized step response gives the fraction of the eventual cascade response observed by time $t$. This tutorial uses that form to specify a likelihood and a practical estimation strategy. For a different known input profile, the corresponding cascade response is obtained by applying the earlier material on convolution responses to that input; the likelihood-fitting principles remain the same.

Definition: Let $F_n(t;k)$ be the normalized step response of a cascade with $n$ identical first-order stages and common rate constant $k>0$. For integer $n\geq 1$, $$ F_n(t;k)=1-\exp(-kt)\sum_{j=0}^{n-1}\frac{(kt)^j}{j!},\qquad t\geq 0. $$ A measurement model for response $y_i$ at time $t_i$ is $y_i=B+G F_n(t_i;k)+\epsilon_i$, where $B$ is baseline, $G$ is the response gain for the sustained input, and $\epsilon_i$ is measurement error.

The likelihood strategy below assumes that the input and sampling times are known, observations are independent conditional on the model, and errors have a common variance and a normal distribution: $\epsilon_i\sim N(0,\sigma^2)$. Thus $y_i\sim N(B+G F_n(t_i;k),\sigma^2)$. The gain and baseline are treated as continuous parameters; the stage number $n$ is a positive integer, not a continuous quantity. The model also assumes the same cascade structure and rate constant apply throughout the observation period. These are modeling assumptions, not guarantees conferred by maximizing a likelihood.

Likelihood and a useful profiling step

For observations $y_1,\ldots,y_m$, the likelihood is the product of their normal densities. Equivalently, maximizing it is the same as minimizing the residual sum of squares for the predicted means, provided the common normal-error assumption holds. For a fixed candidate $n$ and $k$, the values $F_n(t_i;k)$ are known. The prediction is then linear in $B$ and $G$, so these two parameters can be estimated by ordinary least squares for that candidate.

\ell(B,G,k,n,\sigma^2) =-\frac{m}{2}\log(2\pi\sigma^2) -\frac{1}{2\sigma^2}\sum_{i=1}^m \left[y_i-B-GF_n(t_i;k)\right]^2

For fixed $n,k$, let $f_i=F_n(t_i;k)$, and define $\bar f$ and $\bar y$ as the means of the $f_i$ values and observations. When the sampled $f_i$ values are not all equal, the least-squares estimates are

\widehat G(n,k)= \frac{\sum_{i=1}^m(f_i-\bar f)(y_i-\bar y)} {\sum_{i=1}^m(f_i-\bar f)^2}, \qquad \widehat B(n,k)=\bar y-\widehat G(n,k)\bar f.

Otherwise, $B$ and $G$ are not separately identified by these observations.

When the baseline is known to be zero, the fitted gain instead has the simpler through-origin form $\widehat G=\sum_i f_i y_i/\sum_i f_i^2$. In either case, these estimates can be substituted into the residual sum of squares. The remaining nonlinear search is over $k$ and the integer $n$. For each candidate stage number, one can minimize the residual sum of squares over $k>0$, then compare the resulting likelihoods or a penalized criterion such as AIC.

Definition: Under the independent normal-error model, the maximum-likelihood estimate of the variance for a fitted model with $m$ observations is $\widehat{\sigma}^2_{\mathrm{ML}}=\mathrm{SSE}/m$ when $\mathrm{SSE}>0$, where $\mathrm{SSE}=\sum_i(y_i-\widehat y_i)^2$. If $\mathrm{SSE}=0$, the likelihood is unbounded as $\sigma^2\to0$, so there is no finite positive-variance MLE. This likelihood estimate divides by $m$, not by a residual degrees-of-freedom adjustment; an unbiased residual variance estimate, when appropriate, uses a different denominator.

Profiling out $B$ and $G$ reduces the numerical search and avoids asking an optimizer to find all parameters simultaneously when two have simple conditional solutions. It does not make $k$ or $n$ automatically identifiable. A flat likelihood over a range of $k$, or nearly equal fits across stage numbers, signals that the observed time course may not contain enough information to distinguish them.

Worked estimation and stage comparison

Worked Example: Estimating Gain From a Fixed Cascade Shape

Hypothetical data. Suppose baseline is known to be zero, $n=2$, and $k=0.5$ per hour. At five selected times, the normalized response fractions are approximately $f=(0.300000,0.500000,0.593994,0.700000,0.900057)$; the corresponding hypothetical observations, in response units, are $y=(6.000,6.000,6.200,6.000,8.223099)$. The times are chosen to yield those fractions under the stated cascade; the fraction at $t=4$ hours is $1-e^{-2}(1+2)\approx0.593994$.

For a known zero baseline, the gain estimate is $\widehat G=\sum_i f_i y_i/\sum_i f_i^2$. Using the displayed values gives a numerator of approximately $20.084021$ and a denominator of approximately $1.992931$. Therefore, $$ \widehat G=\frac{20.084021}{1.992931}\approx10.0776\ \text{response units}. $$

At $t=4$ hours, the fitted response is $10.0776(0.593994)\approx5.986$ response units. The residual, observed minus fitted, is $6.200-5.986\approx0.214$ response units. This residual is one contribution to the total squared error; it does not by itself establish whether the model is adequate.

Worked Example: Estimating Measurement Error After Fitting

Continue with the same hypothetical observations and fractions. The fitted responses, using $\widehat G=10.0776$, are approximately $3.023,\ 5.039,\ 5.986,\ 7.054,$ and $9.070$ response units. Subtracting these from the five observations gives residuals of approximately $2.977,\ 0.961,\ 0.214,\ -1.054,$ and $-0.847$ response units.

The sum of squared errors is approximately $2.977^2+0.961^2+0.214^2+(-1.054)^2+(-0.847)^2\approx11.66\ (\text{response units})^2$. With $m=5$, the normal-likelihood estimate is $\widehat{\sigma}^2_{\mathrm{ML}}=11.66/5\approx2.332\ (\text{response units})^2$, and $\widehat{\sigma}_{\mathrm{ML}}=\sqrt{2.332}\approx1.527$ response units. This calculation assumes the same independent, normally distributed error variance at all five times. It estimates noise conditional on the specified cascade shape; it cannot show that this shape is correct.

Worked Example: Comparing Two Candidate Stage Numbers

Suppose a separate hypothetical response time course has $m=12$ observations. After profiling baseline and gain and optimizing $k$, fitting a two-stage model gives $\mathrm{SSE}_2=18.0\ (\text{response units})^2$, while a three-stage model gives $\mathrm{SSE}_3=14.0\ (\text{response units})^2$. For illustration, use $\mathrm{AIC}=m\log(\mathrm{SSE}/m)+2p+C$, where $p$ counts estimated parameters and $C$ is the same constant for both models. Each model estimates $B,G,k,\sigma^2$, so $p=4$.

For two stages, $\mathrm{AIC}_2=12\log(18/12)+8+C\approx12\log(1.5)+8+C\approx12.87+C$. For three stages, $\mathrm{AIC}_3=12\log(14/12)+8+C\approx12\log(1.1667)+8+C\approx9.85+C$. The three-stage model has the lower AIC by about $3.02$, so it is favored by this criterion among these two candidates. This is a relative model comparison, not proof that the true cascade has three stages; other candidates, the sampling design, and model assumptions still matter.

Why sampling times determine what can be estimated

The response fraction starts near zero, rises, and approaches one. If all observations are taken very late, then $F_n(t;k)$ is close to one across the sample. The data can then estimate the eventual level $B+G$, but have little leverage to separate gain from baseline or to identify the rate and stage number. If all observations are very early, the response may remain close to zero, making the signal small relative to measurement error. Measurements concentrated only near the plateau can leave the rising phase poorly determined.

A useful schedule therefore spans the expected response: include pre-response or baseline information when baseline is not known, several measurements during the rising portion, observations around the expected maximum slope, and enough later measurements to reveal approach to the sustained-response plateau. Sampling plans should be informed by plausible ranges of $k$ and $n$, not only a single best guess. When feasible, examine predicted time courses across those plausible parameter values before fixing the schedule.

The stage number changes the shape and delay, while $k$ changes the time scale. With incomplete time coverage, these effects can compensate: a model with more stages and a faster rate may resemble one with fewer stages and a slower rate over the observed interval. A practical warning is a broad profile likelihood for $k$, different parameter combinations with nearly identical fitted curves, or stage-number rankings that change markedly when one time point is removed. Those are signs of weak information, even if the optimizer reports precise-looking point estimates.

Fitting workflow, alternatives, and limitations

Specify the response and error model.
State the known input, baseline treatment, observation times, stage-response equation, and assumptions about error distribution and dependence.

Fit each candidate structure.
For each integer $n$ in a scientifically plausible range, search over $k>0$; at each trial value, estimate linear parameters such as $B$ and $G$, then compute the likelihood.

Compare and assess information.
Compare likelihood-based criteria for candidate stage numbers and inspect parameter uncertainty or profile behavior. Check whether measurements cover both the rising and later portions of the response.

A joint numerical optimization over $B,G,k,\sigma$ is an alternative to profiling the linear parameters. It may be convenient for more complex models, but it can be slower and sensitive to starting values. A Bayesian fit is another approach, incorporating prior distributions and producing posterior uncertainty; it still depends on a correct observation model and sufficiently informative data. For a known, non-step input, fitting the appropriate convolution response is preferable to forcing a step-response curve onto the data.

The simple likelihood described here has limitations. Independent, constant-variance normal errors may be implausible for repeated measurements, assay data with proportional error, or observations with a detection limit. Correlated errors or heteroscedasticity require a different likelihood. Likewise, the profile procedure estimates parameters within a chosen cascade family; it cannot make an incorrect input history, baseline specification, or stage architecture correct. Detailed residual assessment for structural misspecification is a separate task.

Common mistake

Common pitfall: “The model with the smallest residual sum of squares must be the correct stage number.” The correction is that adding flexibility generally cannot increase the best-fitting residual sum of squares, so fit alone rewards more complex models. Use a stated comparison rule, report uncertainty, and consider whether the schedule could distinguish the candidates. Even a preferred candidate remains conditional on the assumed input, response structure, and error model.

A related arithmetic pitfall is rounding the response fractions too early or mixing sums from rounded and unrounded values. Compute the likelihood and fitted values from full-precision model predictions, then round the reported estimates. In the gain example, $10.0776\times0.593994\approx5.98603$, not $5.985$; the corresponding residual is approximately $0.214$ response units.

Check Your Understanding

Answer each question, then compare with the solutions.

Why is stage number searched as an integer while the rate constant is optimized as a positive continuous parameter?
In the model with unknown baseline, why are both $B$ and $G$ estimated for each trial $n,k$?
What does the likelihood estimate $\mathrm{SSE}/m$ represent under the stated normal-error model?
Why can late-only sampling make gain and rate difficult to distinguish?
If a three-stage model has a lower AIC than a two-stage model, what conclusion is justified?

Solutions. (1) The number of stages counts discrete stages and cannot be fractional; $k>0$ describes a continuous rate. (2) For fixed $n,k$, the predicted response is linear in baseline and gain, and both affect the fitted level; they can be estimated by least squares before comparing nonlinear parameters. (3) It is the maximum-likelihood estimate of the common error variance, in squared response units, conditional on the fitted model and independent normal errors. (4) Late fractions are all near one, so predictions are nearly $B+G$ and vary little with $k$; the time-course shape is then weakly represented. (5) The three-stage model is preferred by AIC among the compared candidates under the stated criterion, but this does not prove it is the true structure.