  Inherent stochasticity in gene expression leads to distributions of mRNA copy
numbers in a population of identical cells. These distributions are determined
primarily by the multitude of states of a gene promoter, each driving
transcription at a different rate. In an era where single-cell mRNA copy number
data are more and more available, there is an increasing need for fast
computations of mRNA distributions. In this paper, we present a method for
computing separate distributions for each species of mRNA molecules, i. e.
mRNAs that have been either partially or fully processed post-transcription.
The method involves the integration over all possible realizations of promoter
states, which we cast into a set of linear ordinary differential equations of
dimension $M\times n_j$, where $M$ is the number of available promoter states
and $n_j$ is the mRNA copy number of species $j$ up to which one wishes to
compute the probability distribution. This approach is superior to solving the
Master equation (ME) directly in two ways: a) the number of coupled
differential equations in the ME approach is
$M\times\Lambda_1\times\Lambda_2\times ...\times\Lambda_L$, where $\Lambda_j$
is the cutoff for the probability of the $j^{\text{th}}$ species of mRNA; and
b) the ME must be solved up to the cutoffs $\Lambda_j$, which are {\it ad hoc}
and must be selected {\it a priori}. In our approach, the equation for the
probability to observe $n$ mRNAs of any species depends only on the the
probability of observing $n-1$ mRNAs of that species, thus yielding a correct
probability distribution up to an arbitrary $n$. To demonstrate the validity of
our derivations, we compare our results with Gillespie simulations for ten
randomly selected system parameters.
