Stochastic effects in cellular systems are usually modeled and simulated with Gillespie's stochastic simulation algorithm (SSA), which follows the same theoretical derivation as the chemical master equation (CME), but the low efficiency of SSA limits its application to large chemical networks. To improve efficiency of stochastic simulations, Haseltine and Rawlings proposed a hybrid of ODE and SSA algorithm, which combines ordinary differential equations (ODEs) for traditional deterministic models and SSA for stochastic models. In this dissertation, accuracy analysis, efficient implementation strategies, and application of of Haseltine and Rawlings's hybrid method (HR) to a budding yeast cell cycle model are discussed. Accuracy of the hybrid method HR is studied based on a linear chain reaction system, motivated from the modeling practice used for the budding yeast cell cycle control mechanism. Mathematical analysis and numerical results both show that the hybrid method HR is accurate if either numbers of molecules of reactants in fast reactions are above certain thresholds, or rate constants of fast reactions are much larger than rate constants of slow reactions. Our analysis also shows that the hybrid method HR allows for a much greater region in system parameter space than those for the slow scale SSA (ssSSA) and the stochastic quasi steady state assumption (SQSSA) method. Implementation of the hybrid method HR requires a stiff ODE solver for numerical integration and an efficient event-handling strategy for slow reaction firings. In this dissertation, an event-handling strategy is developed based on inverse interpolation. Performances of five wildly used stiff ODE solvers are measured in three numerical experiments. Furthermore, inspired by the strategy of the hybrid method HR, a hybrid of ODE and SSA stochastic models for the budding yeast cell cycle is developed, based on a deterministic model in the literature. Simulation results of this hybrid model match very well with biological experimental data, and this model is the first to do so with these recently available experimental data. This study demonstrates that the hybrid method HR has great potential for stochastic modeling and simulation of large biochemical networks. / Ph. D.

Numerous challenges arise in modeling and simulation as biochemical networks are discovered with increasing complexities and unknown mechanisms. With the improvement in experimental techniques, biologists are able to quantify genes and proteins and their dynamics in a single cell, which calls for quantitative stochastic models for gene and protein networks at cellular levels that match well with the data and account for cellular noise. This dissertation studies a stochastic spatiotemporal model of the Caulobacter crescentus cell cycle. A two-dimensional model based on a Turing mechanism is investigated to illustrate the bipolar localization of the protein PopZ. However, stochastic simulations are often impeded by expensive computational cost for large and complex biochemical networks. The hybrid stochastic simulation algorithm is a combination of differential equations for traditional deterministic models and Gillespie's algorithm (SSA) for stochastic models. The hybrid method can significantly improve the efficiency of stochastic simulations for biochemical networks with multiscale features, which contain both species populations and reaction rates with widely varying magnitude. The populations of some reactant species might be driven negative if they are involved in both deterministic and stochastic systems. This dissertation investigates the negativity problem of the hybrid method, proposes several remedies, and tests them with several models including a realistic biological system. As a key factor that affects the quality of biological models, parameter estimation in stochastic models is challenging because the amount of empirical data must be large enough to obtain statistically valid parameter estimates. To optimize system parameters, a quasi-Newton algorithm for stochastic optimization (QNSTOP) was studied and applied to a stochastic budding yeast cell cycle model by matching multivariate probability distributions between simulated results and empirical data. Furthermore, to reduce model complexity, this dissertation simplifies the fundamental cooperative binding mechanism by a stochastic Hill equation model with optimized system parameters. Considering that many parameter vectors generate similar system dynamics and results, this dissertation proposes a general α-β-γ rule to return an acceptable parameter region of the stochastic Hill equation based on QNSTOP. Different objective functions are explored targeting different features of the empirical data. / Doctor of Philosophy / Modeling and simulation of biochemical networks faces numerous challenges as biochemical networks are discovered with increased complexity and unknown mechanisms. With improvement in experimental techniques, biologists are able to quantify genes and proteins and their dynamics in a single cell, which calls for quantitative stochastic models, or numerical models based on probability distributions, for gene and protein networks at cellular levels that match well with the data and account for randomness. This dissertation studies a stochastic model in space and time of a bacterium’s life cycle— Caulobacter. A two-dimensional model based on a natural pattern mechanism is investigated to illustrate the changes in space and time of a key protein population. However, stochastic simulations are often complicated by the expensive computational cost for large and sophisticated biochemical networks. The hybrid stochastic simulation algorithm is a combination of traditional deterministic models, or analytical models with a single output for a given input, and stochastic models. The hybrid method can significantly improve the efficiency of stochastic simulations for biochemical networks that contain both species populations and reaction rates with widely varying magnitude. The populations of some species may become negative in the simulation under some circumstances. This dissertation investigates negative population estimates from the hybrid method, proposes several remedies, and tests them with several cases including a realistic biological system. As a key factor that affects the quality of biological models, parameter estimation in stochastic models is challenging because the amount of observed data must be large enough to obtain valid results. To optimize system parameters, the quasi-Newton algorithm for stochastic optimization (QNSTOP) was studied and applied to a stochastic (budding) yeast life cycle model by matching different distributions between simulated results and observed data. Furthermore, to reduce model complexity, this dissertation simplifies the fundamental molecular binding mechanism by the stochastic Hill equation model with optimized system parameters. Considering that many parameter vectors generate similar system dynamics and results, this dissertation proposes a general α-β-γ rule to return an acceptable parameter region of the stochastic Hill equation based on QNSTOP. Different optimization strategies are explored targeting different features of the observed data.