Tutorials › Quantitative Systems Pharmacology (QSP) › Deriving Deterministic Limits from Stochastic Cell Models

Delay and Stochastic Extensions · Tutorial 102 of 799

Deriving Deterministic Limits from Stochastic Cell Models

Learn how stochastic transition propensities generate deterministic equations for population densities as population size grows, and when that approximation is informative.

Advanced 10 min read

What You'll Learn

  • Express a population CTMC in terms of count states, transition increments, and propensities.
  • Scale propensities and states to obtain the limiting density drift.
  • State the regularity, initial-condition, and boundedness assumptions behind the deterministic limit.
  • Derive deterministic equations for linear birth–death and logistic population models.
  • Explain why sublinear-rate channels can have zero limiting drift yet matter in finite populations.
  • Diagnose when a deterministic limit may not represent the cell-level behavior of interest.

From Stochastic Transitions to Population-Scale Equations

A cell-population model may specify events one at a time: a cell divides, dies, or changes state at a random time. Such a model describes both the average population trend and its random fluctuations. When the population is large, it is often useful to ask whether the normalized population follows a deterministic rate equation. The answer depends on how transition rates scale with population size, not simply on the fact that the system contains many cells.

The starting point is a continuous-time Markov chain (CTMC) with count state \(X^N(t)\in\mathbb{Z}_{\ge 0}^d\). The index \(N\) denotes a population-size scale. There are finitely many event types, indexed by \(r=1,\ldots,R\). Event \(r\) changes the count vector by a fixed increment \(\nu_r\in\mathbb{Z}^d\), and occurs at propensity \(a_r^N(x)\) when the current count state is \(x\). A propensity is the event rate in the current state, with units of events per unit time. This transition-and-propensity description is the starting point of the earlier material on formulating cellular CTMCs.

Definition: The normalized state, or density, is \(Z^N(t)=X^N(t)/N\). A density-dependent scaling has limiting rate functions \(\alpha_r\) such that, for every compact set of densities, $$ \sup_{z}\left|\frac{a_r^N(Nz)}{N}-\alpha_r(z)\right|\longrightarrow 0 \quad\text{as }N\longrightarrow\infty. $$ The candidate deterministic drift is $$ F(z)=\sum_{r=1}^{R}\nu_r\alpha_r(z), $$ so the corresponding rate equation is \(\frac{dz}{dt}=F(z)\).

The uniform convergence in the definition is understood over feasible states (or their appropriate integer approximations) in the compact set. Each count increment \(\nu_r\) is assumed not to grow with \(N\). Since one event changes the density by \(\nu_r/N\), events occurring at rates of order \(N\) can produce an order-one density change over a fixed time interval. The sum of their expected density changes per unit time is precisely \(F(z)\).

Assumptions Behind the Deterministic Limit

The rate equation is a large-population approximation, not an identity for a finite CTMC. A standard finite-time law of large numbers applies under conditions such as the following:

  • Convergent initial densities: \(Z^N(0)\) converges in probability to a fixed density \(z_0\).
  • Density-scaled rates: the propensities satisfy the uniform-on-compact scaling above, with finitely many event types and fixed increments.
  • Regular limiting drift: \(F\) is locally Lipschitz on the relevant region, so the initial-value problem has a unique solution while that solution remains in the region.
  • Control of escape and explosion: over the time interval of interest, the chains remain in bounded regions with probability approaching one, or suitable growth conditions provide equivalent control. Rates must not produce infinitely many events in finite time.

With these conditions, the density trajectory \(Z^N(t)\) converges in probability, uniformly over each fixed bounded time interval on which the solution remains controlled, to the solution \(z(t)\) of $$ \frac{dz}{dt}=F(z),\qquad z(0)=z_0. $$ The conclusion concerns the density path as \(N\) grows. It does not say that a finite population has no stochastic variation, nor does it automatically give a useful approximation for rare events or for times that grow with \(N\).

Why the Drift Becomes Deterministic

The justification follows from separating the accumulated expected change from the random deviation around it. For the normalized CTMC, each event of type \(r\) contributes a jump \(\nu_r/N\), and its rate at density \(z\) is \(a_r^N(Nz)\). Thus the instantaneous expected density change is $$ \sum_{r=1}^{R}\frac{\nu_r}{N}a_r^N(Nz) =\sum_{r=1}^{R}\nu_r\frac{a_r^N(Nz)}{N} \longrightarrow \sum_{r=1}^{R}\nu_r\alpha_r(z)=F(z). $$

More precisely, the density process can be written as its initial value, an integrated drift, and a martingale \(M^N(t)\):

$$ Z^N(t)=Z^N(0)+\int_0^t F^N(Z^N(s))\,ds+M^N(t), \qquad F^N(z)=\sum_{r=1}^{R}\nu_r\frac{a_r^N(Nz)}{N}. $$

On a bounded region where the propensities are of order \(N\), the rate of accumulated squared density jumps is of order \(N(1/N^2)=1/N\). For example, the predictable quadratic variation of coordinate \(j\) of the martingale over time \(T\) is $$ \left\langle M_j^N\right\rangle(T) =\int_0^T\sum_{r=1}^{R} \frac{(\nu_{r,j})^2}{N^2}a_r^N(X^N(s))\,ds, $$ which is bounded by a constant times \(T/N\) on that region. Therefore the martingale deviation tends to zero in probability on a fixed finite interval. The scaled drift \(F^N\) approaches \(F\); with the stated local regularity, the integral equation then approaches the unique ODE solution. In particular, if the drift error is uniformly at most \(\epsilon_N\) on the region and the Lipschitz constant of \(F\) there is \(L\), the difference between the paths is bounded by the initial discrepancy and martingale discrepancy, plus the accumulated drift error, multiplied by at most \(e^{LT}\). Each term tends to zero under the assumptions, giving the finite-time convergence.

Worked Example: Linear Birth–Death Population

Consider a hypothetical cell population with count \(X^N\). Each cell divides at rate \(\beta=0.6\) per day and dies at rate \(\delta=0.2\) per day. The two increments are \(+1\) and \(-1\), with propensities \(a_+^N(x)=\beta x\) and \(a_-^N(x)=\delta x\). These rates are linear in the count and satisfy the density scaling: $$ \frac{a_+^N(Nz)}{N}=0.6z,\qquad \frac{a_-^N(Nz)}{N}=0.2z. $$

The limiting functions are \(\alpha_+(z)=0.6z\) and \(\alpha_-(z)=0.2z\). Multiplying each event rate by its count increment gives $$ \frac{dz}{dt}=(+1)(0.6z)+(-1)(0.2z)=0.4z \quad\text{per day}. $$

If the initial count is \(400\) in a population scale \(N=1000\), the initial density is \(z_0=0.4\). The deterministic solution is \(z(t)=0.4e^{0.4t}\). At one day, it predicts \(z(1)=0.4e^{0.4}\approx0.5967\), or about \(597\) cells per \(1000\) on this population scale. This is a density prediction; it is not an assertion that a finite CTMC must contain exactly that number of cells.

State-Dependent Rates and Finite-Size Corrections

Density dependence can arise when cells compete, interact, or encounter a limited resource. The large-population equation uses the limiting rate functions; terms that vanish after scaling are not retained in the leading-order drift. It is important to check that scaling directly rather than assume every propensity must be exactly proportional to \(N\).

Worked Example: Logistic Growth from Birth and Competition Death

Suppose births occur at rate \(a_+^N(x)=0.8x\) per day, and deaths occur at rate $$ a_-^N(x)=0.2x+0.6\frac{x(x-1)}{N} \quad\text{per day}. $$ The first death term is a background death rate; the second represents pairwise competition. The factor \(x(x-1)\) counts ordered pairs of distinct cells, while division by \(N\) keeps the total competition rate of order \(N\) when \(x\) is of order \(N\).

At density \(z=x/N\), the birth rate divided by \(N\) is \(0.8z\). The death rate divided by \(N\) is $$ \frac{a_-^N(Nz)}{N} =0.2z+0.6z\left(z-\frac{1}{N}\right) \longrightarrow 0.2z+0.6z^2. $$ Since births add one cell and deaths remove one, the deterministic equation is $$ \frac{dz}{dt}=0.8z-(0.2z+0.6z^2)=0.6z(1-z) \quad\text{per day}. $$

For \(N=1000\) and \(x=500\), the hypothetical birth rate is \(0.8(500)=400\) births per day. The death rate is \(0.2(500)+0.6(500)(499)/1000=100+149.7=249.7\) deaths per day. The net count drift is \(400-249.7=150.3\) cells per day, or \(150.3/1000=0.1503\) density units per day. The limiting drift at \(z=0.5\) is \(0.6(0.5)(1-0.5)=0.15\) per day, close to the finite-\(N\) value.

The comparison also illustrates a diagnostic: finite-size terms may create a small discrepancy between the exact scaled drift and its limit. Here that discrepancy comes from \(1/N\). It vanishes as \(N\) increases. If a term instead remains substantial in the scaled propensity, it belongs in the limiting rate function and must not be discarded.

Worked Example: A Sublinear-Rate Channel

Suppose a modeled event adds one cell and has propensity \(a_{\mathrm{rare}}^N(x)=1\) per hour, independent of the current count while the event is available. Its contribution to the scaled propensity is $$ \frac{a_{\mathrm{rare}}^N(Nz)}{N}=\frac{1}{N}\longrightarrow 0. $$ It therefore satisfies the density-dependent scaling with limiting rate function \(\alpha_{\mathrm{rare}}(z)=0\), and its contribution to the limiting density drift is zero.

That does not make the event impossible or necessarily irrelevant in a finite population. Over two hours, if the event remains available at rate one per hour, the probability of at least one occurrence is \(1-e^{-2}\approx0.8647\). Each occurrence changes density by only \(1/N\), so its direct effect on the density vanishes as \(N\) grows. But an occurrence might still matter to a finite-population question, such as whether a rare cell type is seeded. The deterministic density equation alone does not answer that question.

When to Use the Limit—and When to Be Cautious

A deterministic equation is useful when the target is a population-scale trajectory over a fixed time horizon and the density-scaled assumptions are plausible. It offers a compact description of the leading-order trend and can be easier to analyze than a distribution of stochastic paths. It is not a replacement for the CTMC when finite-count events, variability, or the probability of a rare outcome is central.

ApproachWhat it representsMain limitation
Deterministic rate equationLimiting density trajectory as \(N\) growsDoes not retain finite-population randomness or zero-limit-rate events
Exact CTMC simulationIndividual stochastic event histories under the specified CTMCOne trajectory is variable; many runs may be needed for distributional summaries
Tau-leapingAn approximate finite-population simulation using event-count stepsIntroduces time-step approximation error and requires appropriate step sizes

The CTMC and deterministic model should be compared on the quantity of interest. Agreement of typical density trajectories does not establish agreement in extinction probability, time to a rare event, or the frequency of a low-count phenotype. Nor does a good fit at one population scale establish that the assumed scaling is correct at other scales.

Useful checks include evaluating \(a_r^N(Nz)/N\) over relevant densities and several values of \(N\), checking whether it approaches a stable function, and confirming that the proposed ODE remains in a biologically feasible region. If scaled propensities fail to converge, or if the limiting drift is discontinuous or non-unique in a relevant region, the standard deterministic-limit argument may not apply as stated. Sublinear-rate channels are not a failure merely because their limit is zero; a channel fails the scaling condition only when its scaled propensity does not converge to the required limiting function.

Common Mistake

Common pitfall: “Every channel must have a propensity exactly proportional to \(N\), so a constant-rate event is excluded.” This is incorrect. A propensity equal to \(1\) has scaled value \(1/N\to0\), so it satisfies the scaling with a zero limiting rate contribution. The correction is to test convergence of \(a_r^N(Nz)/N\). A zero limit means that the channel does not contribute to the leading-order density drift; it does not mean the finite-population event cannot occur or matter.

A second related mistake is to interpret the ODE as the exact expected count trajectory. The deterministic-limit result concerns convergence of normalized paths under its assumptions. It does not, by itself, prove that the expected CTMC state equals the ODE solution at finite \(N\).

Check Your Understanding

Answer each question, then compare with the solutions below.

  1. A channel has increment \(-1\) and propensity \(a^N(x)=0.3x\) per day. What is its limiting contribution to the density drift?
  2. A channel has increment \(+1\) and propensity \(a^N(x)=2\) per day. Does it satisfy the scaling, and what is its limiting drift contribution?
  3. For \(N=500\), a population has \(x=100\), a birth propensity \(0.5x\) per day, and a death propensity \(0.1x\) per day. Calculate the current density and the scaled net drift, then give the limiting drift at that density.
  4. Name two conditions needed to interpret the rate equation as a finite-time large-population limit.

Solutions. (1) The limiting rate is \(\alpha(z)=0.3z\). Multiplying by the increment gives \((-1)(0.3z)=-0.3z\) density units per day. (2) The scaled propensity is \(2/N\to0\), so the channel satisfies the scaling with \(\alpha(z)=0\), and its limiting density drift contribution is zero. (3) The density is \(100/500=0.2\). The count drift is \(0.5(100)-0.1(100)=40\) cells per day, so the scaled net drift is \(40/500=0.08\) density units per day. The limiting expression is \((0.5-0.1)z=0.4(0.2)=0.08\) per day. (4) Examples include convergence in probability of the initial densities, uniform-on-compact convergence of scaled propensities, a locally Lipschitz limiting drift, and control of escape or explosion on the time interval considered.

Clinical Biostats