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 speci cally 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.