With the sequence of genomes of many organisms now available, the major chal- lenge of functional genomics is 'reassembling the pieces'. A chemical reaction network is considered to be a very simple and ecient view of a living system. The goal here is to be able to simulate an arbitrary ensemble of hypothesized biochemical and gene regulatory networks to predict what a cell is doing and compare them with the observed data. Any biological network can be modeled as a system of Ordinary Dif- ferential Equations (ODEs) mathematically. There are a number of standard ODE solvers, such as the Euler and Runge Kutta (RK) methods. But from time to time, these methods could be very inecient for some special systems, called sti sys- tems. We consider other special ODE solvers, such as Livermore Solver for ODE with General Sparse Jacobian Matrices (LSODES). A simulator KINSOLVER with 4 integration options (Euler, Modied Euler, RK, Adaptive RK-Fehlberg) is now updated by adding LSODES. I also presented the theoretical and numerical work on the stiness of some biological circuits like qa gene cluster, lac operon and oregonator. Performance of dierent methods including MATLAB ODE suite are compared too. The eciency of computation could be improved by a factor of 10 compared to the standard 4th order RK method.