accepted by J. Comput. Biology 07/24/2009Cells of all organisms share the ability to respond to various extracellular signals. Depending on the cell type and the organism, these signals may include hormones secreted by other cells or changes in nutrient concentrations. The signals are processed by an intricate network of protein-protein interactions, including phosphorylation and de-phosphorylation events. As some signaling proteins are only present in low concentrations, random fluctuations may affect the dynamics of the network. The mathematical modeling of networks with significant random fluctuations requiresthe use of stochastic methods. The stochastic dynamics of a chemical reaction system are described by the Chemical Master Equation. Often the numerical evaluation of this equation is problematic. The first problem is that many systems have an infinite number of possible states; leaving simulations of individualtrajectories as the only alternative. To circumvent this problem, we focus on a class of systems that have a finite state space. More specifically, we focus on networks of phosphorylation cycles without taking into account the synthesis and degradation of proteins. Thesecond problem is that memory requirements cause a practical limit to the size of systems that can be evaluated. In this paper, we discuss how these limitations can be overcome using parallel computation and methods dealing efficiently with the availablememory. These methods were implemented in a parallel C++ program. We discuss two networks for which the stochastic dynamics were evaluated using this program: a single phosphorylation cycle and an oscillating MAP-kinase cascade.