In this thesis we study stochastic modeling and simulation methods for biochemical systems. The thesis is focused on systems with multi-state and multi-scale features and divided into two parts. In the first part, we propose new algorithms that improve existing multi-state simulation methods. We first compare the well known Gillespie\'s stochastic simulation algorithm (SSA) with the StochSim, an agent-based simulation method. Based on the analysis, we propose a hybrid method that possesses the advantages of both methods. Then we propose two new methods that extend the Network-Free Algorithm (NFA) for rule-based models. Numerical results are provided to show the performance improvement by our new methods. In the second part, we investigate two simulation schemes for the multi-scale feature: Haseltine and Rawlings\' hybrid method and the quasi-steady-state stochastic simulation method. We first propose an efficient partitioning strategy for the hybrid method and an efficient way of building stochastic cell cycle models with this new partitioning strategy. Then, to understand conditions where the two simulation methods can be applied, we develop a way to estimate the relaxation time of the fast sub-network, and compare it with the firing interval of the slow sub-network. Our analysis are verified by numerical experiments on different realistic biochemical models.