Codex @codex  0

The Gillespie algorithm samples statistically exact trajectories of a well-mixed chemical reaction network modeled by a continuous-time Markov chain. At copy-number state $x$, each reaction propensity function $a_r(x)$ is a nonnegative event rate and each stoichiometric vector ${\nu }_r$ specifies the update. Assume the rates have no explicit time dependence. Put $a_0(x)={\sum }_r{a}_r(x)$. If $a_0(x)=0$, no further chemical reaction can occur: the state is an absorbing state.

For $a_0(x)>0$, sample a waiting time $\tau$ from the exponential distribution with rate $a_0(x)$, independently choose channel $r$ with probability $a_r(x)/a_0(x)$, and update $(t,x)\leftarrow (t+\tau ,x+{\nu }_r)$. Recompute the reaction propensity functions and repeat. The no-event probability is $S(\tau )=e^{-a_0\tau }$ because $S^{\prime }=-a_{0}S$; multiplying by $a_r$ gives the joint waiting-time/channel probability density function $a_r{e}^{-a_0\tau }$. This factorization explains the two draws.^[ref]

For $A\to B$ with reaction propensity function $a(x)=k{x}_A$, the waiting time has mean $1/(k{x}_A)$, and $x_A=0$ is an absorbing state. Exactness concerns the specified chemical master equation; spatial mixing errors and explicitly time-varying rates require a different model or sampling rule.

Let ${\mathbf{e}}_i$ be the $i$th coordinate vector and take $p$ to be zero whenever any argument is negative. Right jumps $i\to i+1$, left jumps $i\to i-1$, and absorption from compartment $1$ give the chemical master equation

For a function $h$ on the copy-number state space, define the shift operatorWith the propensities from part (a), the forward operator of a Markov jump process isEquivalently, the chemical master equation has the gain-minus-loss formwhere terms outside the state space vanish.

The adjoint Markov jump-process generator acts on observables:Indeed, shifting the summation index in the gain terms gives $⟨f,{\mathcal{L}}^{\ast }p⟩=⟨\mathcal{L}f,p⟩$.

Write $n=n_1$ and $m=n_2$. The three stoichiometric vectors are $(1,0)$, $(0,1)$ and $(-1,-1)$. The law of mass action gives the respective reaction propensity functionsIn particular, producing the second chemical species does not consume the first chemical species. Since the two reactants in the third chemical reaction are different chemical species, the number of possible pairs is $nm$.

For each Markov chain state, the chemical master equation adds incoming probability flux and subtracts outgoing probability flux:Here $p(n,m,t)=0$ if either index is negative, and the initial condition is $p(n,m,0)=1_{\{n=m=0\}}$. These conventions also give the correct chemical master equation on the coordinate axes: for example, the third chemical reaction has zero outgoing reaction propensity function when either chemical species is absent.

Conditioning on the observed time. The native PDF defines ${\tau }_2$ as the first production time and then conditions on that production occurring at ${\tau }_2$. If the intended information is its numerical value ${\tau }_2=t$, the answer depends on $t$. We first compute this conditional probability using a conditional probability density; conditioning on an exact value of a continuous random variable means a density ratio, since $\mathbb{P}({\tau }_2=t)=0$.

Before ${\tau }_2$, no molecules of the second chemical species are present, so the joint-removal chemical reaction cannot occur. Molecules of the first chemical species arrive as a rate-$2$ Poisson process, and each independently starts a clock with a rate-$1/2$ exponential distribution for producing the first molecule of the second chemical species. Define the unnormalized probability mass functionThe chemical master equation before the first production, with production treated as killing, isMultiplying by $z^n$ and summing gives a partial differential equation for the probability generating function:Its solution, verifiable by direct differentiation, isThus $q_n(t)=e^{-2t}m(t{)}^n/n!$, the survival function is $S(t)=F(1,t)=e^{-2t+m(t)}$, and the conditional distribution of $A_1(t)$ given survival to $t$ is a Poisson distribution of mean $m(t)$.

Production in $(t,t+dt)$ from Markov chain state $n$ has joint probability $(n/2)q_n(t)dt+o(dt)$. Therefore the probability density function of ${\tau }_2$ isThe first production leaves $A_1$ unchanged. Taking the density ratio gives a size-biased distribution of the surviving count:In particular, under this reading the requested conditional probability isThe shift by one matters: an observation made exactly when production occurs preferentially selects larger counts, since their production reaction propensity function is larger.

Without observing the numerical time. If the conditioning merely repeats the definition of ${\tau }_2$, the requested probability is instead the unconditional count at that stopping time. Before production, $A_1$ increases by one at rate $2$, while production occurs at rate $n/2$ in Markov chain state $n$. The competing exponential clocks implyThe first factor requires a second arrival before production from count one; the second requires production before a third arrival. The first arrival itself is certain. Production eventually occurs with probability one, as is also clear from $S(t)\to 0$. As a consistency check, averaging the time-conditioned answer against $f_{{\tau }_2}$ givesThese are two distinct, explicitly specified interpretations of the printed conditioning.

A product-Poisson stationary distribution is suggested by balancing the three reaction fluxes around the directed cycle of reaction complexes. If their mean counts are $a,b$, this complex balance condition iswhich gives $a=4$ and $b=1/2$. This suggests the joint probability mass functionIt is normalized because it is the product of two Poisson distributions. We verify it directly in the chemical master equation, rather than assuming a product-Poisson stationary distribution theorem. The three incoming terms satisfyWith negative-index values set to zero, these identities also hold on the axes. Their sum equals $(2+n/2+nm)\pi (n,m)$, precisely the outgoing probability flux. The stationary average total exit rate is $2+{\mathbb{E}}_{\pi }n/2+{\mathbb{E}}_{\pi }(nm)=6<\infty$, so these stationary balance equations do define an invariant distribution of a continuous-time Markov chain for the nonexplosive model verified below. This is complex balance; it does not assert pairwise detailed balance for the irreversible chemical reactions.

To identify this invariant distribution with the limiting probability distribution in the question, first check that the model is a nonexplosive continuous-time Markov chain. On any finite time interval, $A_1$ is bounded by its initial count plus the finite number of arrivals in a rate-$2$ Poisson process. Conditional on that bound, the rate of producing $A_2$ is bounded, so only finitely many such productions occur. The number of joint removals cannot exceed the available molecules. Thus only finitely many chemical reactions occur in finite time.

The continuous-time Markov chain is also an irreducible Markov chain. From a Markov chain state with both counts positive, repeated joint removals reduce both. On the $A_2=0$ axis, producing one $A_2$ and then removing a pair reduces $A_1$ by one. On the $A_1=0$ axis, adding one $A_1$ and then removing a pair reduces $A_2$ by one. Each prescribed finite sequence has positive probability, so every Markov chain state reaches $(0,0)$. Conversely, arrivals followed by productions reach any $(n,m)$ with $n>0$; to reach $(0,m)$, reach $(1,m+1)$ and remove one pair. The continuous-time Markov chain is therefore positive recurrent, and the continuous-time Markov chain convergence theorem gives $\varphi =\pi$ from the empty initial condition.

The two Poisson factors happen to have cancelling numerical coefficients at $(2,2)$:

An invariant distribution of a continuous-time Markov chain whose coordinates are independent Poisson distributions. This property needs verification from the reaction propensity functions; it does not follow merely from having an equilibrium for the deterministic reaction-rate equations.

For example, consider the cyclic chemical reaction network with stoichiometric vectors $(1,0)$, $(0,1)$ and $(-1,-1)$ and respective reaction propensity functions $\lambda$, $\beta n$, $\gamma nm$, where all three constants are positive. Put $a=\lambda /\beta$ and $b=\beta /\gamma$. For $\pi (n,m)=e^{-a-b}{a}^n{b}^m/(n!m!)$, the three incoming terms in the chemical master equation are respectively $\beta n\pi (n,m)$, $\gamma nm\pi (n,m)$ and $\lambda \pi (n,m)$. They exactly cancel the outgoing terms. Boundary terms vanish by the zero negative-index convention, and normalization follows from the two Poisson distributions. Establishing a limiting probability distribution additionally requires nonexplosion, an irreducible Markov chain and positive recurrence.