  In this work, we consider the problem of estimating summary statistics to
characterise biochemical reaction networks of interest. Such networks are often
described using the framework of the Chemical Master Equation (CME). For
physically-realistic models, the CME is widely considered to be analytically
intractable. A variety of Monte Carlo algorithms have therefore been developed
to explore the dynamics of such networks empirically. Amongst them is the
multi-level method, which uses estimates from multiple ensembles of sample
paths of different accuracies to estimate a summary statistic of interest. {In
this work, we develop the multi-level method in two directions: (1) to increase
the robustness, reliability and performance of the multi-level method, we
implement an improved variance reduction method for generating the sample paths
of each ensemble; and (2) to improve computational performance, we demonstrate
the successful use of a different mechanism for choosing which ensembles should
be included in the multi-level algorithm.
