  While the analytical solution for the marginal distribution of a stochastic
chemical reaction network has been extensively studied, its joint distribution,
i.e. the solution of a high-dimensional chemical master equation, has received
much less attention. Here we develop a novel method of computing the exact
joint distributions of a wide class of first-order stochastic reaction systems
in steady-state conditions. The effectiveness of our method is validated by
applying it to four gene expression models of biological significance,
including models with 2A peptides, nascent mRNA, gene regulation, translational
bursting, and alternative splicing.
