Biochemical systems have important practical applications, in particular to understanding key intra-cellular processes. Often biochemical kinetic models represent cellular processes as systems of chemical reactions, traditionally modelled by the deterministic reaction rate equations. In the cellular environment, many biological processes are inherently stochastic. The stochastic fluctuations due to the presence of some low molecular populations may have a great impact on the biochemical system behavior. Then, stochastic models are required for an accurate description of the system dynamics. Biochemically reacting systems often evolve on multiple time-scales, thus their mathematical models manifest stiffness. Stochastic models which, in addition, are stiff are computationally challenging, therefore the need for developing effective and accurate numerical methods for approximating their solution. We present an adaptive stepsize strategy for the numerical solution of stochastic models of biochemical systems. Numerical experiments on models of practical interest show the advantages of the adaptive strategies over the fixed-step methods.