Why the Source of Variation Matters
Repeated simulations of a cellular model can produce different outcomes for two distinct reasons. Even with the same parameter values and initial conditions, discrete reaction events can generate different trajectories. This is demographic noise. Outcomes can also differ because the model parameters differ, for example across plausible parameter sets. That is parameter-driven variation in the model prediction.
These sources answer different questions. Variation among stochastic trajectories at fixed parameters describes the spread generated by discrete events under that specified model. Variation across parameter values describes how the model’s expected outcome changes when the parameter vector changes. The latter depends on which parameter values are considered plausible or relevant; it is not automatically a property of the cell population.
The target outcome must be specified before either component is summarized. It could be the cell count at a fixed time, the time to extinction, or an indicator that a threshold was crossed. A variance for one outcome and time point does not describe all trajectory variability.
The Variance Decomposition
Suppose parameter vectors are drawn from a specified distribution, and conditional on each vector the model generates a stochastic outcome. The overall variance across both sources of variation satisfies the law of total variance:
The first term averages the conditional demographic variance over the parameter distribution. The second term measures how much the conditional mean outcome changes across that distribution. This is a decomposition of the variance for the specified mixture of parameter values and stochastic trajectories; changing the parameter distribution can change the second term and the total.
A short derivation makes clear why the components add. Write \(\mu(\theta)=\operatorname{E}[Y\mid\theta]\) and \(\mu=\operatorname{E}[Y]\). Then \(Y-\mu=(Y-\mu(\theta))+(\mu(\theta)-\mu)\). Squaring and taking expectations gives two squared terms and a cross-term. Conditional on \(\theta\), \(\operatorname{E}[Y-\mu(\theta)\mid\theta]=0\), so the cross-term has expectation zero. Therefore:
This identity does not require the conditional variance to be constant across parameter values. It does require a well-defined outcome with finite variance under the joint sampling scheme. In an application, the interpretation also requires stating how parameter vectors are selected. A distribution representing uncertainty about a fixed but unknown parameter has a different meaning from a distribution representing genuine differences among cells or patients.
A Nested Simulation Design
A direct design samples parameter vectors in an outer loop and simulates multiple stochastic trajectories at each vector in an inner loop. Let \(Y_{ij}\) be the outcome from trajectory \(j\) at parameter vector \(\theta_i\), with \(i=1,\ldots,I\) and \(j=1,\ldots,M\). Use the same initial conditions and outcome definition throughout unless their variation is deliberately part of the question.
For each parameter vector, calculate the mean \(\bar{Y}_i\) and sample variance \(s_i^2\) across its \(M\) trajectories. The average of the \(s_i^2\) values estimates average conditional demographic variance for the sampled parameter vectors. The sample variance of the \(\bar{Y}_i\) values describes the spread in conditional means across those vectors, but each \(\bar{Y}_i\) is itself estimated with simulation error. Replication within each parameter vector helps separate that error from genuine differences in conditional means.
Worked Example: Estimating Both Components From Nested Simulations
Suppose three hypothetical parameter vectors are each simulated for three independent stochastic trajectories. The outcome is the number of viable cells at 24 hours.
| Parameter vector | Three simulated outcomes (cells) | Mean (cells) | Sample variance (cells²) |
|---|---|---|---|
| A | 8, 10, 12 | 10 | 4 |
| B | 15, 17, 19 | 17 | 4 |
| C | 22, 24, 26 | 24 | 4 |
Within each parameter vector, the deviations from its mean are \(-2,0,2\) cells. Thus, for A, the sample variance is \((4+0+4)/(3-1)=4\) cells²; the same calculation applies to B and C. Their average within-vector sample variance is \(4\) cells², an estimate of average demographic variance over these three vectors.
The conditional means are 10, 17, and 24 cells, with grand mean \((10+17+24)/3=17\) cells. Their sample variance is:
Thus the estimated between-vector spread in conditional means is 49 cells², while the estimated average within-vector spread is 4 cells². The parameter-driven component is much larger in this small illustration. These are estimates based on three parameter vectors and three trajectories per vector, not exact population components; in particular, the between-vector estimate is affected by simulation error in each mean.
Pairing Simulations to Compare Parameter Settings
Pairing is useful when the aim is to compare outcomes under two or more parameter settings. For each replicate index, simulate each setting using a valid shared source of stochastic randomness, then analyze the within-pair differences. If paired outputs are positively correlated, common stochastic fluctuations can partly cancel in the difference. This can make the contrast between settings easier to estimate than a contrast based on unrelated trajectories.
The pairing does not remove demographic noise from either setting, and it does not change either setting’s marginal outcome distribution. It changes the covariance between simulated outputs and therefore the variability of their difference. The random-number construction must preserve the correct marginal stochastic process for each setting. Simply restarting a pseudorandom generator at the same seed is not always sufficient: if event sequences diverge, the apparent pairing may be arbitrary or misleading.
Worked Example: A Paired Contrast Across Two Parameter Settings
Consider a hypothetical comparison of two parameter settings for the same 24-hour viable-cell outcome. Four paired simulations use a coupling designed to preserve both settings’ marginal stochastic distributions. The outcomes, in cells, are:
| Pair | Setting A | Setting B | Difference, B − A |
|---|---|---|---|
| 1 | 8 | 14 | 6 |
| 2 | 10 | 17 | 7 |
| 3 | 12 | 20 | 8 |
| 4 | 9 | 16 | 7 |
The mean paired difference is \(\bar{d}=(6+7+8+7)/4=7\) cells. The sample standard deviation of the differences is:
The standard error of the mean paired difference, as a Monte Carlo summary of these four paired replicates, is \(s_d/\sqrt{4}\approx0.4082\) cells. The estimated contrast is 7 cells, with setting B higher in these simulated pairs. This calculation does not establish that the settings differ by 7 cells in a biological population; it summarizes a small simulation comparison. Pairing is helpful here because the paired differences vary less than the individual outcomes, but that benefit depends on the coupling and the correlation it induces.
When Demographic Variance Depends on Parameters
Parameter changes can affect both the conditional mean and the conditional variance. For example, a parameter setting that increases the number of reaction events may also increase trajectory-to-trajectory dispersion. The law of total variance still applies: average the within-setting variances using the parameter distribution, then add the variance of the conditional means. Do not assume a single common demographic variance unless the model or diagnostics support that simplification.
Worked Example: Unequal Demographic Variance Across Settings
In another hypothetical nested simulation, three parameter vectors give endpoint cell counts of 2, 4, 6; 8, 12, 16; and 18, 20, 22, respectively. The within-vector means are 4, 12, and 20 cells. The sample variances are \(4\), \(16\), and \(4\) cells²: for the middle vector, for example, the squared deviations sum to \(16+0+16=32\), and \(32/(3-1)=16\) cells².
The average estimated demographic variance is \((4+16+4)/3=8\) cells². The sample variance of the three conditional means is:
The two estimated components are therefore 8 and 64 cells². The noticeably larger within-vector variance at the middle setting is a warning against describing demographic noise as constant across parameters. The decomposition remains interpretable, but the average within-vector component conceals this variation; report the individual conditional variances as well when that pattern matters.
Checks, Alternatives, and Limitations
Before interpreting a decomposition, check that the design matches the question. The outcome, time point, initial conditions, and parameter distribution should be explicit. At each parameter vector there must be enough independent stochastic trajectories to estimate conditional spread. The parameter vectors must cover the distribution for which the between-vector component is intended. If some settings have much larger conditional variance, inspect those settings rather than relying only on an average.
A common mistake is to call the variance of all simulated outcomes “demographic noise.” That pooled variance combines within-parameter stochastic variation with variation in conditional means across the parameter values used. The correction is to preserve parameter identity in the simulation output, estimate within-vector and between-vector components separately, and state the parameter distribution behind the latter.
A deterministic simulation at each parameter vector can show how conditional mean predictions vary with parameters, but it cannot represent demographic variance. Conversely, many stochastic trajectories at one fixed parameter vector estimate demographic variation there, but cannot characterize parameter-driven variation. A nested design addresses both. Paired simulations are an efficient comparison tool for contrasts when valid common randomness induces useful correlation; they are not a substitute for estimating the variance components in the chosen joint model.
Finally, parameter uncertainty is not automatically biological heterogeneity. If parameters are uncertain because evidence is limited, the second component describes uncertainty in predictions under the chosen parameter distribution. If parameters genuinely differ across cells, it may describe population heterogeneity. Correlations among parameter values, misspecified parameter distributions, finite simulation counts, and model structural error can all change the conclusions. Neither component captures uncertainty about whether the model structure itself is adequate.
Check Your Understanding
Answer each question using the stated outcome and parameter design. Brief solutions follow.
- At one fixed parameter vector, 20 stochastic trajectories produce an endpoint count distribution. Which variance component can these simulations estimate?
- Why does pooling trajectories from several parameter vectors and calculating one sample variance fail to isolate demographic noise?
- For three equally weighted parameter settings, the conditional means are 5, 8, and 11 units. Calculate their sample variance.
- What must be true of a shared-randomness construction before paired simulation differences can be interpreted?
- Can parameter-driven variation in predicted means be interpreted as biological heterogeneity whenever a parameter distribution is used? Explain.
Solutions. (1) They estimate demographic variance conditional on that parameter vector, but not variation across parameter values. (2) The pooled variance includes both within-vector trajectory variation and between-vector differences in conditional means. (3) The mean is 8 units, so the sample variance is \(((5-8)^2+(8-8)^2+(11-8)^2)/(3-1)=9\) units². (4) It must induce dependence for variance reduction while preserving the correct marginal stochastic distribution for each setting. (5) No. If the distribution represents uncertainty about a fixed unknown parameter, the component describes prediction uncertainty; biological heterogeneity is appropriate only when the parameter variation represents genuine differences among biological units.