Reaction networks are usually modeled by continuous time Markov chains. It is well known that in regimes where many reactions happen in short time intervals, exact path simulation methods, like the SSA, may be computationally infeasible. On the other hand, the approximate tau-leap method, based on the random time change representation by Kurtz, has the undesirable feature of leading to negative (non-physical) population numbers with positive probability. In this talk, we present a Chernoff type bound for setting an upper bound for this exit probability and controlling the step size of the tau-leap. We propose a hybrid exact/approximate adaptive method for path simulation that avoids unnecessarily small tau-leap step sizes. Next, we show how to couple two hybrid paths associated to two consecutive levels in a hierarchy of time discretizations. This is the key point to define a Multilevel Monte Carlo estimator of the expected value of an observable of the process. We provide the simulation setting that approximately optimizes the computational cost while controlling the accuracy of the method.