Why first-passage times matter
Two stochastic cell-population models can have similar average cell counts over time but very different probabilities of rapid extinction. In some applications, the clinically or biologically relevant question is not just how many cells are expected at a particular time, but when a population first reaches a specified state. This timing can be variable from one realization of the same model to another.
A first-passage distribution describes that variation across simulated trajectories. It is not the trajectory of an average population, and a threshold crossing of the average trajectory is not a substitute for it. The earlier material on moment equations concerns summaries such as means; those summaries alone generally do not determine the distribution of a hitting time.
Define the target before measuring its time
Let \(X(t)\) denote the stochastic state at time \(t\). For a cell-count model, \(X(t)\) may be a count or a vector of counts. A target set \(A\) specifies which states count as having reached the event of interest. The first-passage time is the earliest time the process enters that set.
The target and the population being counted must be stated precisely. For example, extinction of one cell subtype means that subtype's count first reaches zero; it does not necessarily mean that every cell in a multi-type system is gone. If the initial state is already in \(A\), then \(\tau_A=0\). A threshold target can be crossed without being absorbing: a population can later grow above \(k\), but its first-passage time remains the time of its first entry into \(N(t)\leq k\).
For a fixed model, parameter set, and initial condition, let \(\tau_A\) vary across independent simulated trajectories. Its cumulative distribution function (CDF) is \(F_A(t)=P(\tau_A\leq t)\), and its survival function is \(S_A(t)=P(\tau_A>t)=1-F_A(t)\). In this context, “survival” means not yet having reached the target; it need not mean that the cells are biologically healthy. If some trajectories never reach the target, the distribution can have probability mass at \(\infty\), and \(F_A(t)\) need not approach 1 as \(t\) increases.
Worked Example: Reading Extinction and Threshold Times From One Path
Suppose a hypothetical exact stochastic simulation gives the following cell counts. The times shown are event times at which the count changes.
| Time (hours) | Cell count after the event |
|---|---|
| 0.0 | 12 |
| 0.7 | 10 |
| 1.1 | 7 |
| 1.6 | 4 |
| 2.4 | 2 |
| 3.0 | 0 |
The path is piecewise constant between events. For extinction, the target is \(A=\{0\}\), so this trajectory has \(\tau_A=3.0\) hours. For threshold \(k=5\), the first state with \(N(t)\leq5\) occurs at 1.6 hours, giving \(\tau_5=1.6\) hours. For threshold \(k=3\), the first such state occurs at 2.4 hours, so \(\tau_3=2.4\) hours. These are three different event definitions applied to the same trajectory.
With an event-driven simulation, a jump that lands in the target identifies the crossing time at that event. If results are recorded only at a coarse set of observation times, the exact crossing may not be known; it is only known to lie between observations. Reporting the later observation time as an exact event time would overstate the timing precision.
Estimate the distribution from complete trajectories
To estimate a first-passage distribution, simulate multiple independent trajectories under the same model, parameters, initial condition, and target definition. For trajectory \(i\), record its first-passage time \(\tau_i\). When every trajectory reaches the target within the observation window, the empirical CDF at time \(t\) is the fraction of recorded times no greater than \(t\):
Here \(n\) is the number of simulated trajectories and \(I(\tau_i\leq t)\) equals 1 when trajectory \(i\) has reached the target by time \(t\), and 0 otherwise. The empirical CDF is a step function. It estimates the probability of reaching the specified target by each time, under the model and simulation conditions used. It is not a confidence statement about whether a particular trajectory will cross.
Worked Example: Empirical Extinction Distribution
Consider eight hypothetical independent trajectories, all starting from the same cell count. Their recorded extinction times, in hours, are \(2,3,3,5,6,6,8,10\). All eight reached zero before the simulation ended, so the complete-data empirical CDF applies.
By 4 hours, three of the eight trajectories have become extinct. Thus \(\widehat{F}_A(4)=3/8=0.375\). By 5 hours, four have become extinct: \(\widehat{F}_A(5)=4/8=0.500\). The estimated probability of extinction by 5 hours is therefore 0.500 for this simulated model setup.
The empirical 25th percentile, using the definition \(\inf\{t:\widehat{F}_A(t)\geq0.25\}\), is 3 hours: the CDF first reaches 0.25 at that time. The empirical median under the same lower-quantile definition is 5 hours because the CDF first reaches 0.50 at 5 hours. The sample mean extinction time is \((2+3+3+5+6+6+8+10)/8=43/8=5.375\) hours. The mean and median summarize different aspects of the same simulated times; neither replaces the full stepwise distribution.
Trajectories that have not reached the target
A path that has not reached the target when simulation stops is right-censored: its first-passage time is known to exceed its observed follow-up time, but its eventual time is unknown. If a trajectory runs to a common fixed horizon \(T\) without reaching the target, its record is “no event by \(T\),” not an event time equal to \(T\) and not evidence that \(\tau_A=\infty\). Treating such records as if the event occurred at \(T\) biases the estimated distribution toward earlier events.
For right-censored data, the Kaplan–Meier estimator uses the risk set at each distinct event time. The risk set contains trajectories still under observation and not yet known to have reached the target immediately before that time. At an event time \(t_j\), let \(n_j\) be the number at risk and \(d_j\) the number of events. The estimated survival probability is the product of the conditional fractions that remain event-free at each event time.
Censored trajectories reduce later risk sets but do not count as events. The estimator requires event times to be observed accurately and censoring to be independent of the event time, at least conditional on any relevant groups or covariates. Administrative censoring at a fixed simulation horizon meets this condition when all trajectories are followed to that horizon and the horizon is not determined by their event risk. A state-dependent stopping rule can violate it. The estimate describes the distribution only over times with adequate follow-up; it cannot reveal what happens beyond the observed range.
Worked Example: Kaplan–Meier Estimate With Censored Paths
Five hypothetical trajectories are followed until extinction or their respective stopping times. The records are: extinction at 2 hours; no extinction by 3 hours; extinction at 4 hours; no extinction by 5 hours; extinction at 6 hours. Assume censoring is independent of extinction time.
At 2 hours, all five trajectories are at risk and one becomes extinct. The estimated survival fraction is multiplied by \(1-1/5=4/5=0.800\). The trajectory censored at 3 hours is then removed from the later risk sets; it is not counted as an extinction.
At 4 hours, three trajectories remain at risk and one becomes extinct. The estimated survival is now \(0.800(1-1/3)=0.800(2/3)=0.5333\). The trajectory censored at 5 hours is removed next. At 6 hours, one trajectory remains at risk and becomes extinct, so the product becomes \(0.5333(1-1/1)=0\).
Thus the Kaplan–Meier estimates of extinction by 2 and 4 hours are \(1-0.800=0.200\) and \(1-0.5333=0.4667\), respectively. The estimated median is 6 hours, the first event time at which estimated extinction probability reaches at least 0.5. The final estimate of zero survival is a consequence of the last observed at-risk trajectory having an event; with so few paths, it is not strong evidence that every trajectory in the underlying model eventually becomes extinct.
Justification and interpretation
The empirical CDF follows directly from the event indicator. For a complete trajectory, \(I(\tau_i\leq t)\) is 1 exactly when the trajectory has reached the target by \(t\). Averaging these indicators therefore gives the observed fraction that has reached it. For censored records, the indicator at times after censoring is unknown, so simply averaging observed event indicators would incorrectly treat unknown outcomes as no event.
The Kaplan–Meier product addresses this by estimating conditional survival at each event time. Among the \(n_j\) trajectories still at risk just before \(t_j\), \(d_j\) events occur, leaving an estimated conditional survival fraction \(1-d_j/n_j\). Multiplying these conditional fractions gives the estimated probability of remaining event-free through all event times up to \(t\). Subtracting from 1 gives the estimated cumulative probability of passage by \(t\). This construction depends on censoring not selectively removing trajectories with different unobserved event risks.
Define which cells are counted and the target set, such as total count equal to zero or a specified count at or below \(k\).
For every simulated trajectory, record the first-passage time if observed; otherwise record its follow-up time and mark it as right-censored.
Use the empirical CDF for complete event times or the Kaplan–Meier estimate for right-censored times. State the horizon, units, event definition, and number of trajectories.
Diagnostics, alternatives, and limitations
Check that the recorded time is the first entry, not a later visit to the target. Confirm that the target definition is consistent across trajectories and that all simulations use the intended initial condition and parameter set. If a simulation reports only regularly spaced observations, assess whether a crossing may have occurred between recorded times. In leap-based simulations, a step can jump across a threshold; the step endpoint may not identify the exact first-passage time. Report the time resolution or use a simulation method that resolves the event timing when exact timing matters.
The complete-data empirical CDF is simple and transparent when every trajectory reaches the target. The Kaplan–Meier estimator is the standard alternative when some trajectories are right-censored under an independent-censoring assumption. Reporting only a mean or median is more compact, but it conceals distributional shape and tail behavior. A deterministic trajectory or a crossing time computed from the mean cell count answers a different question: it describes the model’s mean state, not the distribution of stochastic first-passage times.
Finite simulations estimate, rather than exactly determine, the model-implied distribution. Estimates at times beyond the longest meaningful follow-up are unsupported, and an apparent absence of events by the horizon does not establish that the event is impossible. Likewise, a point estimate of an extreme quantile can be unstable when few trajectories inform that part of the distribution. These estimates describe stochastic variability conditional on the model and its fixed inputs; they do not by themselves account for parameter uncertainty or establish that the model represents a biological system accurately.
Check Your Understanding
Answer each question using the event definition and observation information given.
- A path starts with 20 cells and first reaches 4 cells at 2.5 hours, then later grows to 8. What is its first-passage time to \(N(t)\leq5\)?
- Six complete extinction times are 1, 2, 2, 4, 5, and 8 hours. What is the empirical extinction probability by 2 hours?
- A simulation stops at 10 hours with a positive cell count. Should its extinction time be recorded as 10 hours, infinity, or right-censored at 10 hours?
- Why can the crossing time of a mean cell-count trajectory not generally be interpreted as the median extinction time?
Solutions. (1) The first state at or below 5 occurs at 2.5 hours, so the first-passage time is 2.5 hours; later growth does not change it. (2) Three of the six times are no greater than 2 hours, so \(\widehat{F}_A(2)=3/6=0.500\). (3) It is right-censored at 10 hours: the event has not been observed, and its later occurrence or non-occurrence is unknown. (4) The mean trajectory averages cell counts at each time, whereas the median extinction time is a quantile of trajectory-level first-passage times. Averaging paths does not preserve their crossing-time distribution.