Constructing Reaction-Network CTMC Rates

From Reaction Channels to State Changes

A reaction network specifies possible events, but a continuous-time Markov chain (CTMC) records states. To pass from one description to the other, identify how each reaction changes the state, calculate its propensity at the current state, and combine rates when different reactions produce the same new state. The resulting rates determine both the generator and the holding-time and jump-destination laws.

Let \((\Omega,\mathcal{F},\mathbb{P})\) be a probability space, and let \((X_t)_{t\ge0}\) be a continuous-time process with values in a state space \(S\subseteq\mathbb{Z}_{\ge0}^{d}\). Its natural filtration is \(\mathcal{F}_t=\sigma(X_s:0\le s\le t)\). A state \(x=(x_1,\ldots,x_d)\) records the counts of \(d\) species. A reaction channel \(r\) has a reactant vector \(\alpha_r\), a product vector \(\beta_r\), and net change \(\nu_r=\beta_r-\alpha_r\). If the channel fires at state \(x\), the new state is \(x+\nu_r\).

The state space must contain the states under consideration and be closed under every feasible state-changing reaction. In this tutorial, assume the total rate of state-changing reactions at each state being considered is finite. A propensity gives a rate for a reaction channel; it is not automatically the rate of a particular state-to-state move, because channels can have the same net change.

Constructing the Generator

Why rates are added. Suppose two distinct channels both move the process from \(x\) to \(y\). In a short interval of length \(h\), either channel can produce that move. The probability that channel \(r\) fires is \(a_r(x)h+o(h)\); the probability of two or more firings is \(o(h)\) when the total rate is finite. Adding the single-channel probabilities gives \((a_1(x)+a_2(x))h+o(h)\), so the state-to-state rate is the sum, not either channel’s rate alone.

This reasoning also verifies the generator row sum. Every state-changing channel at \(x\) belongs to exactly one destination \(y\ne x\), so summing \(q_{xy}\) over those destinations counts each such channel’s propensity exactly once. Therefore \(\sum_{y\ne x}q_{xy}=q_x\), and adding \(q_{xx}=-q_x\) gives zero.

Holding Times and Jump Destinations

At a fixed state \(x\), the reaction-channel construction assigns an exponential clock with rate \(a_r(x)\) to each channel. Channels with zero rate cannot fire. The first clock to ring determines the next reaction. If multiple channels cause the same state change, the state process records their common destination, not which channel fired.

Proof. Assume first \(q_x>0\). For each positive-rate channel \(r\), let \(E_r\) be its exponential clock. Independence of the clocks gives, for \(t\ge0\), \(\mathbb{P}(H>t\mid X_0=x)=\prod_r\mathbb{P}(E_r>t)=\prod_r e^{-a_r(x)t}=e^{-q_xt}\). Zero-rate channels may be omitted because they never ring. This is the survival function of an exponential random variable with rate \(q_x\).

For a positive-rate channel \(r\), the probability density that its clock rings at time \(t\) while every other clock exceeds \(t\) is \(a_r(x)e^{-q_xt}\). Consequently, for a destination \(y\ne x\), adding the densities of exactly the channels with \(x+\nu_r=y\) gives the joint density \(q_{xy}e^{-q_xt}\) for “the first state change occurs at time \(t\) and goes to \(y\).” Integrating over \(t\ge0\) gives \(\mathbb{P}(X_H=y\mid X_0=x)=q_{xy}\int_0^\infty e^{-q_xt}\,dt=q_{xy}/q_x\). The joint density factors as \(q_{xy}e^{-q_xt}= (q_{xy}/q_x)(q_xe^{-q_xt})\), the product of the destination probability and the density of \(H\). This factorization proves independence. If \(q_x=0\), every state-changing rate is zero, so no state-changing clock rings; the state stays at \(x\). \(\square\)

In particular, the jump probabilities are obtained by dividing a row’s off-diagonal rates by its exit rate. They describe the destination conditional on a jump, not the speed of the process: multiplying all rates in a row by the same positive constant leaves those probabilities unchanged but changes the holding-time distribution.

Interpreting the Construction Correctly

A generator is a state-based summary of reaction propensities. It keeps the total rate of each possible state change, but it does not necessarily preserve the identity of the reaction channel. If channel labels are themselves part of the model’s observable output, they can be recorded separately; the CTMC on species counts alone only records the resulting state.

A reaction with zero net change, such as a channel whose product and reactant vectors are equal, does not produce a jump of the count-valued state process. Its firing may be meaningful in a richer model that records reaction events, but it is not a transition to a different state and is therefore excluded from \(q_x\). In contrast, an absorbing state is a state at which every feasible state-changing propensity is zero; there the state process has no future jumps.

Common mistake. In the second example, it is wrong to assign the rate from \(k=0\) to \(k=1\) as \(1\), using only the channel \(A\to B\), and then set the diagonal to \(-1\). The other channel also takes \(k=0\) to \(k=1\), at rate \(2\). The correct rate is \(1+2=3\), so the correct generator row is \((-3,3,0)\), not \((-1,1,0)\). The incorrect row would predict an exponential holding time with rate \(1\), rather than the rate \(3\) implied by both available channels.

For a practical construction, list the channels and their propensities at each state of interest, discard channels that are infeasible or do not change the tracked state, group the remaining rates by destination, and set each diagonal to the negative sum of its off-diagonal entries. Then check that every row sums to zero. The competing-clock theorem supplies the corresponding holding-time and jump-destination laws.