Formulating Cellular CTMCs with State Transitions and Propensities

From Cellular Events to a Stochastic Model

A cellular model may describe events such as a molecule binding a receptor, a cell dividing, or a protein being degraded. In a deterministic model, rates govern continuous changes in amounts or concentrations. A continuous-time Markov chain (CTMC) instead represents the system’s current state as counts and each event as a discrete change to those counts. The event rate depends on the state at that time.

Formulating a CTMC requires more than writing down a list of biological events. For every event, the model needs a precise change in state and a propensity: the instantaneous rate at which that event occurs, conditional on the current state. These definitions determine the stochastic model before any simulation algorithm is chosen.

Definition: A cellular CTMC has a state \(X(t)=(X_1(t),\ldots,X_d(t))\), where each component is a nonnegative integer count. Event channel \(j\) has a transition vector \(\nu_j\), which gives the change in state if that event occurs, and a propensity \(a_j(x)\), which gives its instantaneous event rate when the current state is \(x\). The event changes \(x\) to \(x+\nu_j\).

The state should contain enough information to determine the future event rates. In particular, the Markov assumption is that, conditional on the present state, the distribution of future states does not depend on the earlier path. This assumption is appropriate when the listed counts and any other included state variables fully summarize what matters for future events. If event rates depend on an omitted quantity, such as time since a molecule entered a compartment, the chosen state may not satisfy this assumption.

States, Transition Vectors, and Propensities

Suppose there are \(d\) tracked species or cell classes and \(m\) event channels. The state space is a subset of \(\mathbb{Z}_{\geq 0}^d\). For each channel \(j\), the transition vector \(\nu_j\) has one entry per state component. A negative entry denotes consumption, a positive entry denotes production, and zero denotes no change in that component. The propensity must be nonnegative and must be zero wherever the event would lead to an invalid state.

A propensity has units of inverse time, such as events per minute. Its defining short-interval interpretation is that, for a sufficiently small interval \(h\), the probability of one event of type \(j\) is approximately \(a_j(x)h\), given that the state at the start of the interval is \(x\). More precisely, for a time-homogeneous CTMC with finite rates at \(x\),

P(\text{one event of type } j \text{ in } [t,t+h]\mid X(t)=x) = a_j(x)h+o(h).

Here \(o(h)\) denotes a term whose ratio to \(h\) tends to zero as \(h\) tends to zero. The probability of two or more events in that interval is also \(o(h)\). Thus propensities are rates, not probabilities; multiplying a propensity by a sufficiently short time interval gives a first-order approximation to an event probability.

For a reaction network, event channels are often written as reactions, with the transition vector giving the net count change. For example, \(A+B\rightarrow C\) has changes \(-1\) in \(A\), \(-1\) in \(B\), and \(+1\) in \(C\). The propensity specifies how frequently that event occurs at each possible state. Different channels may lead to the same destination state; their rates then contribute to the total rate of that state change.

Worked Example: First-Order Conversion

Consider conversion of an unbound molecule \(A\) into \(A^*\). Let the state be \(x=(A,A^*)\), and suppose each \(A\) converts independently at rate \(k=0.30\) per minute. Specify the event channel and calculate its propensity at \(x=(4,1)\).

One conversion consumes one \(A\) and produces one \(A^*\), so the transition vector is \(\nu=(-1,+1)\). The propensity is the per-molecule rate multiplied by the number of molecules available:

a(x)=kA=(0.30\ \text{min}^{-1})(4)=1.20\ \text{min}^{-1}.

If this event occurs, the state becomes \(x+\nu=(3,2)\). If \(A=0\), then \(a(x)=0\), so the channel cannot produce an invalid negative count. The propensity of \(1.20\ \text{min}^{-1}\) is an instantaneous rate; for example, over a sufficiently short interval of \(0.001\) minutes, the first-order probability of one conversion is approximately \(1.20(0.001)=0.00120\).

Worked Example: Binding Between Two Species

Suppose \(A\) and \(B\) bind to form complex \(C\), with event \(A+B\rightarrow C\). Use the state order \(x=(A,B,C)\), and define the propensity as \(a(x)=kAB\), where \(k=0.020\) molecule\(^{-1}\) minute\(^{-1}\). At \(x=(5,3,1)\), find the transition and rate.

The event consumes one molecule of each free species and produces one complex, so \(\nu=(-1,-1,+1)\). Substitution gives

a(x)=(0.020\ \text{molecule}^{-1}\text{min}^{-1})(5\ \text{molecules})(3\ \text{molecules}) =0.300\ \text{min}^{-1}.

If the event occurs, the state becomes \((4,2,2)\). The units cancel to give a rate in inverse minutes. If either free-species count is zero, the propensity is zero, as required for this binding event. This example uses a stated convention for \(k\); a model must define its rate-constant convention rather than assume that all reported bimolecular constants use the same units.

Worked Example: Dimerization and Pair Counting

Consider dimer formation \(2A\rightarrow D\), with state \(x=(A,D)\). Suppose the propensity is defined as \(a(x)=c\binom{A}{2}\), where \(c=0.050\) per minute per available pair. At \(x=(4,0)\), calculate the rate and the state after one event.

There are \(\binom{4}{2}=4(3)/2=6\) distinct pairs of \(A\) molecules. The transition consumes two \(A\) molecules and produces one \(D\), so \(\nu=(-2,+1)\). Therefore,

a(x)=(0.050\ \text{min}^{-1}\text{pair}^{-1})(6\ \text{pairs}) =0.300\ \text{min}^{-1}, \qquad x+\nu=(2,1).

The pair-counting form guarantees that the propensity is zero for \(A=0\) or \(A=1\). It also makes the rate-constant convention explicit. Writing a different formula, such as a constant times \(A(A-1)\), is possible only if the constant is scaled consistently with that formula.

Worked Example: Birth and Death of Cells

Let \(N\) be the number of cells in a population. Suppose each cell divides at rate \(\beta=0.40\) per day and dies at rate \(\delta=0.10\) per day. At \(N=6\), specify both channels and their rates.

Division adds one cell, so its transition is \(+1\) and its propensity is \(a_{\mathrm{div}}(N)=\beta N\). Death removes one cell, so its transition is \(-1\) and its propensity is \(a_{\mathrm{death}}(N)=\delta N\). At \(N=6\),

a_{\mathrm{div}}(6)=(0.40\ \text{day}^{-1})(6)=2.40\ \text{day}^{-1}, \qquad a_{\mathrm{death}}(6)=(0.10\ \text{day}^{-1})(6)=0.600\ \text{day}^{-1}.

The total event rate is \(2.40+0.600=3.00\ \text{day}^{-1}\). A division changes the state from 6 to 7; a death changes it from 6 to 5. At \(N=0\), both propensities are zero, so neither channel can leave the allowed nonnegative state space.

Generator and Forward Equation

The transition rates define the local dynamics of the CTMC. For a function \(f\) of the state, the generator describes the instantaneous rate of change in its expected value, conditional on the current state. The expression follows directly from the short-time event probabilities.

Definition: For event channels \(j=1,\ldots,m\), the generator applied to \(f\) is

\mathcal{L}f(x)=\sum_{j=1}^{m}a_j(x)\bigl[f(x+\nu_j)-f(x)\bigr].

To justify this formula, condition on the state \(X(t)=x\). In an interval \(h\), channel \(j\) occurs with probability \(a_j(x)h+o(h)\); with no event, the state remains \(x\). Since two or more events have probability \(o(h)\), the conditional expected change in \(f\) is

E[f(X(t+h))-f(X(t))\mid X(t)=x] =h\sum_{j=1}^{m}a_j(x)\bigl[f(x+\nu_j)-f(x)\bigr]+o(h).

Dividing by \(h\) and taking the limit as \(h\) tends to zero gives the generator formula. This derivation assumes finite total propensity at \(x\), with the event channels and rates well-defined there.

The same rates give a probability balance equation. Let \(p(x,t)=P(X(t)=x)\). Probability flows into \(x\) from states that can reach it in one event and flows out of \(x\) through events originating there:

\frac{d}{dt}p(x,t) = \sum_{j=1}^{m} a_j(x-\nu_j)p(x-\nu_j,t) - \sum_{j=1}^{m} a_j(x)p(x,t).

In the incoming sum, terms with an invalid predecessor \(x-\nu_j\) are treated as zero. This equation is the forward equation (or chemical master equation) for the specified channels. It makes clear why both the transition vectors and propensities are necessary: the vectors identify possible sources and destinations, while the propensities determine the probability flow.

Assumptions, Diagnostics, and Alternatives

A CTMC formulation is appropriate when the modeled events are discrete, their rates are meaningful at the chosen time scale, and the state contains the information needed to determine those rates. For a time-homogeneous model, the propensities depend on the current state and fixed parameters, not explicitly on time. If rates vary with time, the model can instead be specified as a time-inhomogeneous CTMC, with propensities written \(a_j(x,t)\). That change should be explicit.

Before using a formulation, check it against the following conditions:

State validity: every allowed event leaves all counts nonnegative. Rates for impossible events must be zero at the relevant boundary states.
Event accounting: each transition vector matches the event’s net change in every tracked component. A molecule used in an event should not disappear from the bookkeeping unless the model intentionally omits its destination.
Rate units: every propensity has units of inverse time. The rate constant’s units must be compatible with the count factors in its propensity.
Markov sufficiency: the state must include variables needed to determine future rates. Dependence on unrecorded history can invalidate the CTMC state description.
Finite rates: the total propensity must be finite at each state considered. If rates grow without bound as the process evolves, additional analysis may be needed to establish that the model does not make infinitely many transitions in finite time.

For an instantaneous event model, distinct channels are assumed not to occur simultaneously. This is expressed by the short-time property that the probability of two or more events is \(o(h)\). Delays or other history-dependent mechanisms require care: a delayed event cannot generally be represented by an instantaneous transition whose rate depends only on the current counts. One may need additional state variables or a different model formulation.

A deterministic ordinary differential equation is a useful alternative when the aim is to describe continuous changes in population amounts or concentrations. A CTMC instead retains integer counts and random event-to-event variability. The two formulations may use related biological mechanisms, but their state variables and interpretations differ. Transit-chain models, discussed in the earlier material on distributed delays, represent stages through additional compartments; a CTMC can also use discrete counts for such stages, provided each transition and propensity is specified. Neither representation removes the need to check whether the state captures the relevant history.

Common pitfall: Treating a propensity as the probability of an event, or using a deterministic rate expression without checking how it applies to integer counts, can produce an invalid model. A propensity is an instantaneous rate in inverse-time units. Verify its short-interval interpretation, its boundary behavior, and the state change it accompanies.

A good formulation is auditable: a reader can start from any allowed state, identify each possible event, compute its propensity, and verify the resulting state. If that cannot be done unambiguously, the model is not yet specified well enough for simulation or analysis.

Check Your Understanding

Answer each question using the state, transition, and propensity definitions above.

For \(A\rightarrow B\), with state order \((A,B)\), what is the transition vector?
If one molecule of \(A\) converts at rate \(0.2\) per minute and \(A=7\), what is the propensity?
For \(2A\rightarrow D\), why must the propensity be zero when \(A=1\)?
What is wrong with calling a propensity of \(0.5\) per minute a 50% probability of an event?
If an event rate depends on how long a molecule has already spent in a compartment, what must be checked about the state representation?

Solutions. (1) \(\nu=(-1,+1)\), because the event consumes one \(A\) and produces one \(B\). (2) \(a=0.2(7)=1.4\ \text{min}^{-1}\). (3) One molecule cannot supply the two molecules required by the event, so no dimerization transition is possible and the propensity must be zero. (4) The propensity is a rate, not a probability; for a short interval \(h\) minutes, the first-order event probability is approximately \(0.5h\), provided \(h\) is sufficiently small. (5) Check whether the current state includes enough information to determine the rate. If time already spent in the compartment matters but is omitted, the Markov assumption for that state is not satisfied.