  Biochemical reaction networks are often modelled using discrete-state,
continuous-time Markov chains. System statistics of these Markov chains usually
cannot be calculated analytically and therefore estimates must be generated via
simulation techniques. There is a well documented class of simulation
techniques known as exact stochastic simulation algorithms, an example of which
is Gillespie's direct method. These algorithms often come with high
computational costs, therefore approximate stochastic simulation algorithms
such as the tau-leap method are used. However, in order to minimise the bias in
the estimates generated using them, a relatively small value of tau is needed,
rendering the computational costs comparable to Gillespie's direct method.
  The multi-level Monte Carlo method (Anderson and Higham, Multiscale Model.
Simul. 10:146-179, 2012) provides a reduction in computational costs whilst
minimising or even eliminating the bias in the estimates of system statistics.
This is achieved by first crudely approximating required statistics with many
sample paths of low accuracy. Then correction terms are added until a required
level of accuracy is reached. Recent literature has primarily focussed on
implementing the multi-level method efficiently to estimate a single system
statistic. However, it is clearly also of interest to be able to approximate
entire probability distributions of species counts. We present two novel
methods that combine known techniques for distribution reconstruction with the
multi-level method. We demonstrate the potential of our methods using a number
of examples.
