Chemical master equation

atropy computes the dynamics of reaction networks in the stochastic description. Let us first describe what a chemical reaction network is. A reaction network consists of a closed vessel with volume $V$ and chemical species \(S_0, \dots, S_{d-1}\). The population number of a species \(S_i\) is denoted as \(x_i \in \mathbb{N}_0\), where \(x=(x_0, \dots, x_{d-1}) \in \Omega\) is called the state. The state of mixture is changed by \(M\) chemical reactions \(R_0, \dots, R_{M-1}\). Such a reaction \(R_{\mu}\) is given by stoichiometric coefficients \(s_{\mu, i}, r_{\mu, i}\in \mathbb{N}_0\) and has the form

\[s_{\mu, 0} S_0 + \dots + s_{\mu, d-1} S_{d-1} \rightarrow r_{\mu, 0} S_0 + \dots + r_{\mu, d-1} S_{d-1}.\]

The change of state $x$ under such a reaction $R_{\mu}$ is then given by $x+\nu_{\mu}$, where $\nu_{\mu} =r_{\mu} - s_{\mu}$ is the stoichiometric vector and $r_{\mu} = (r_{\mu,0}, \dots, r_{\mu,d-1})$, $s_{\mu} = (s_{\mu,0}, \dots, s_{\mu,d-1})$. There are two different settings which specify how such a state evolves in time. We will consider both of them briefly.

In the deterministic description we need the population number to be differentiable and time dependent. Each reaction $R_{\mu}$ has a reaction rate, which allows us to form rate equations for the species $S_i$ on the product side of the corresponding reaction. If we put the rate equations for the $M$ different reactions together, we end up with a coupled system of ODEs. Even though we can now obtain concentrations for all future times, this approach has a number of limitations. The population number takes real values, which makes no physical sense. We can rewrite our equations using the average of our population numbers, but this, in some situations, leads to a result different from the real average. This is due to the assumption that we can neglect stochastic fluctuations. Thus we want to consider the second approach.

For the stochastic description we assume that the integral population numbers $x \in \Omega \subseteq \mathbb{N}_0^d$ are a homogeneous Markov process. Our aim is now to predict the probability of finding a specific amount of species at some given time. In the previous approach we instead wanted to find the precise concentration at future time points. Using the stochastic description we can formulate the kinetic Chemical master equation

\[\partial_t P(t,x) = \sum\limits_{\mu = 0}^{M-1} (\alpha_{\mu}(x-\nu_{\mu})P(t,x-\nu_{\mu})-\alpha_{\mu}(x)P(t,x)),\]

where the population change of a given reaction $R_{\mu}$ is denoted by the stoichiometric vector $\nu_{\mu}$. We call $\alpha_{\mu}(x)$ the propensity function. $P(t,x)$ is a probability distribution. Note that the kinetic CME is a gain-loss equation, where the first term describes a gain and the second term a loss.