Numerical analysis of stochastic biochemical reaction networks

Numerical solution of the chemical master equation for stochastic reaction networks typically suffers from the state space explosion problem due to the curse of dimensionality and from stiffness due to multiple time scales. The dimension of the state space equals the number of molecular species involved in the reaction network and the size of the system of differential equations equals the number of states in the corresponding continuous-time Markov chain, which is usually enormously huge and often even infinite. Thus, efficient numerical solution approaches must be able to handle huge, possibly infinite and stiff systems of differential equations efficiently. In this thesis, we present efficient techniques for the numerical analysis of the biochemical reaction networks. We present an approximate numerical integration approach that combines a dynamical state space truncation procedure with efficient numerical integration schemes for systems of ordinary differential equations including adaptive step size selection based on local error estimates. We combine our dynamical state space truncation with the method of conditional moments, and present the implementation details and numerical results. We also incorporate ideas from importance sampling simulations into a non-simulative numerical method that approximates transient rare event probabilities based on a dynamical truncation of the state space. Finally, we present a maximum likelihood method for the estimation of the model parameters given noisy time series measurements of molecular counts. All approaches presented in this thesis are implemented as part of the tool STAR, which allows to model and simulate the biochemical reaction networks. The efficiency and accuracy is demonstrated by numerical examples.Numerische Lösungen der chemischen Master-Gleichung für stochastische Reaktionsnetzwerke leiden typischerweise an dem Zustandsraumexplosionsproblem aufgrund der hohen Dimensionalität und der Steifigkeit durch mehrfache Zeitskalen. Die Dimension des Zustandsraumes entspricht der Anzahl der molekularen Spezies von dem Reaktionsnetzwerk und die Größe des Systems von Differentialgleichungen entspricht der Anzahl der Zustände in der entsprechenden kontinuierlichen Markov-Kette, die in der Regel enorm gross und oft sogar unendlich gross ist. Daher müssen numerische Methoden in der Lage sein, riesige, eventuell unendlich grosse und steife Systeme von Differentialgleichungen effizient lösen zu können. In dieser Arbeit beschreiben wir effiziente Methoden für die numerische Analyse biochemischer Reaktionsnetzwerke. Wir betrachten einen inexakten numerischen Integrationsansatz, bei dem eine dynamische Zustandsraumbeschneidung und ein Verfahren mit einem effizienten numerischen Integrationsschema für Systeme von gewöhnlichen Differentialgleichungen benutzt werden. Wir kombinieren unsere dynamische Zustandsraumbeschneidungsmethode mit der Methode der bedingten Momente und beschreiben die Implementierungdetails und numerischen Ergebnisse. Wir benutzen auch Ideen des importance sampling für eine nicht-simulative numerische Methode, die basierend auf der Zustandsraumbeschneidung die Wahrscheinlichkeiten von seltenen Ereignissen berechnen kann. Schließlich beschreiben wir eine Maximum-Likelihood-Methode für die Schätzung der Modellparameter bei verrauschten Zeitreihenmessungen von molekularen Anzahlen. Alle in dieser Arbeit beschriebenen Ansätze sind in dem Software-Tool STAR implementiert, das erlaubt, biochemische Reaktionsnetzwerke zu modellieren und zu simulieren. Die Effizienz und die Genauigkeit werden durch numerische Beispiele gezeigt

Approximation Techniques for Stochastic Analysis of Biological Systems

There has been an increasing demand for formal methods in the design process of safety-critical synthetic genetic circuits. Probabilistic model checking techniques have demonstrated significant potential in analyzing the intrinsic probabilistic behaviors of complex genetic circuit designs. However, its inability to scale limits its applicability in practice. This chapter addresses the scalability problem by presenting a state-space approximation method to remove unlikely states resulting in a reduced, finite state representation of the infinite-state continuous-time Markov chain that is amenable to probabilistic model checking. The proposed method is evaluated on a design of a genetic toggle switch. Comparisons with another state-of-the-art tool demonstrate both accuracy and efficiency of the presented method

Noise-Induced Min Phenotypes in E. coli

The spatiotemporal oscillations of the Escherichia coli proteins MinD and MinE direct cell division to the region between the chromosomes. Several quantitative models of the Min system have been suggested before, but no one of them accounts for the behavior of all documented mutant phenotypes. We analyzed the stochastic reaction-diffusion kinetics of the Min proteins for several E. coli mutants and compared the results to the corresponding deterministic mean-field description. We found that wild-type (wt) and filamentous (ftsZ( −)) cells are well characterized by the mean-field model, but that a stochastic model is necessary to account for several of the characteristics of the spherical (rodA(−)) and phospathedylethanolamide-deficient (PE(−)) phenotypes. For spherical cells, the mean-field model is bistable, and the system can get trapped in a non-oscillatory state. However, when the intrinsic noise is considered, only the experimentally observed oscillatory behavior remains. The stochastic model also reproduces the change in oscillation directions observed in the spherical phenotype and the occasional gliding of the MinD region along the inner membrane. For the PE(−) mutant, the stochastic model explains the appearance of randomly localized and dense MinD clusters as a nucleation phenomenon, in which the stochastic kinetics at low copy number causes local discharges of the high MinD(ATP) to MinD(ADP) potential. We find that a simple five-reaction model of the Min system can explain all documented Min phenotypes, if stochastic kinetics and three-dimensional diffusion are accounted for. Our results emphasize that local copy number fluctuation may result in phenotypic differences although the total number of molecules of the relevant species is high

Approximation Techniques for Stochastic Analysis of Biological Systems

There has been an increasing demand for formal methods in the design process of safety-critical synthetic genetic circuits. Probabilistic model checking techniques have demonstrated significant potential in analyzing the intrinsic probabilistic behaviors of complex genetic circuit designs. However, its inability to scale limits its applicability in practice. This chapter addresses the scalability problem by presenting a state-space approximation method to remove unlikely states resulting in a reduced, finite state representation of the infinite-state continuous-time Markov chain that is amenable to probabilistic model checking. The proposed method is evaluated on a design of a genetic toggle switch. Comparisons with another state-of-art tool demonstrates both accuracy and efficiency of the presented method

Partial differential equations for self-organization in cellular and developmental biology

Understanding the mechanisms governing and regulating the emergence of structure and heterogeneity within cellular systems, such as the developing embryo, represents a multiscale challenge typifying current integrative biology research, namely, explaining the macroscale behaviour of a system from microscale dynamics. This review will focus upon modelling how cell-based dynamics orchestrate the emergence of higher level structure. After surveying representative biological examples and the models used to describe them, we will assess how developments at the scale of molecular biology have impacted on current theoretical frameworks, and the new modelling opportunities that are emerging as a result. We shall restrict our survey of mathematical approaches to partial differential equations and the tools required for their analysis. We will discuss the gap between the modelling abstraction and biological reality, the challenges this presents and highlight some open problems in the field

Novel domain expansion methods to improve the computational efficiency of the Chemical Master Equation solution for large biological networks

Background: Numerical solutions of the chemical master equation (CME) are important for understanding the stochasticity of biochemical systems. However, solving CMEs is a formidable task. This task is complicated due to the nonlinear nature of the reactions and the size of the networks which result in different realizations. Most importantly, the exponential growth of the size of the state-space, with respect to the number of different species in the system makes this a challenging assignment. When the biochemical system has a large number of variables, the CME solution becomes intractable. We introduce the intelligent state projection (ISP) method to use in the stochastic analysis of these systems. For any biochemical reaction network, it is important to capture more than one moment: this allows one to describe the system’s dynamic behaviour. ISP is based on a state-space search and the data structure standards of artificial intelligence (AI). It can be used to explore and update the states of a biochemical system. To support the expansion in ISP, we also develop a Bayesian likelihood node projection (BLNP) function to predict the likelihood of the states. Results: To demonstrate the acceptability and effectiveness of our method, we apply the ISP method to several biological models discussed in prior literature. The results of our computational experiments reveal that the ISP method is effective both in terms of the speed and accuracy of the expansion, and the accuracy of the solution. This method also provides a better understanding of the state-space of the system in terms of blueprint patterns. Conclusions: The ISP is the de-novo method which addresses both accuracy and performance problems for CME solutions. It systematically expands the projection space based on predefined inputs. This ensures accuracy in the approximation and an exact analytical solution for the time of interest. The ISP was more effective both in predicting the behavior of the state-space of the system and in performance management, which is a vital step towards modeling large biochemical systems

Efficient Sampling in Stochastic Biological Models

Even when the underlying dynamics are known, studying the emergent behavior of stochastic biological systems in silico can be computationally intractable, due to the difficulty of comprehensively sampling these models. This thesis presents the study of two techniques for efficiently sampling models of complex biological systems. First, the weighted ensemble enhanced sampling technique is adapted for use in sampling chemical kinetics simulations, as well as spatially resolved stochastic reaction-diffusion kinetics. The technique is shown to scale to large, cell-scale simulations, and to accelerate the sampling of observables by orders of magnitude in some cases. Second, I study the free energy estimates of peptides and proteins using Markov random fields. These graphical models are constructed from physics-based forcefields, uniformly sampled at different densities in dihedral angle space, and free energy estimates are computed using loopy belief propagation. The effect of sample density on the free energy estimates provided by loopy belief propagation is assessed, and it is found that in most cases a modest increase in sample density leads to significant improvement in convergence. Additionally, the approximate free energies from loopy belief propagation are compared to statistically exact computations and are confirmed to be both accurate and orders of magnitude faster than traditional methods in the models assessed