  The stochastic dynamics of biochemical networks are usually modelled with the
chemical master equation (CME). The stationary distributions of CMEs are seldom
solvable analytically, and numerical methods typically produce estimates with
uncontrolled errors. Here, we introduce mathematical programming approaches
that yield approximations of these distributions with computable error bounds
which enable the verification of their accuracy. First, we use semidefinite
programming to compute increasingly tighter upper and lower bounds on the
moments of the stationary distributions for networks with rational
propensities. Second, we use these moment bounds to formulate linear programs
that yield convergent upper and lower bounds on the stationary distributions
themselves, their marginals and stationary averages. The bounds obtained also
provide a computational test for the uniqueness of the distribution. In the
unique case, the bounds form an approximation of the stationary distribution
with a computable bound on its error. In the non-unique case, our approach
yields converging approximations of the ergodic distributions. We illustrate
our methodology through several biochemical examples taken from the literature:
Schl\"ogl's model for a chemical bifurcation, a two-dimensional toggle switch,
a model for bursty gene expression, and a dimerisation model with multiple
stationary distributions.
