  Stochastic modeling of reaction networks is a framework used to describe the
time evolution of many natural and artificial systems, including, biochemical
reactive systems at the molecular level, viral kinetics, the spread of epidemic
diseases, and wireless communication networks, among many other examples. In
this work, we present a novel multilevel Monte Carlo method for kinetic
simulation of stochastic reaction networks that is specifically designed for
systems in which the set of reaction channels can be adaptively partitioned
into two subsets characterized by either "high" or "low" activity. Adaptive in
this context means that the partition evolves in time according to the states
visited by the stochastic paths of the system. To estimate expected values of
observables of the system at a prescribed final time, our method bounds the
global computational error to be below a prescribed tolerance, $TOL$, within a
given confidence level. This is achieved with a computational complexity of
order $O(TOL^{-2})$, the same as with an exact method, but with a smaller
constant. We also present a novel control variate technique based on the
stochastic time change representation by Kurtz, which may dramatically reduce
the variance of the coarsest level at a negligible computational cost. Our
numerical examples show substantial gains with respect to the standard
Stochastic Simulation Algorithm (SSA) by Gillespie and also our previous hybrid
Chernoff tau-leap method.
