Solving the chemical master equation using sliding windows

Background  

The chemical master equation (CME) is a system of ordinary differential equations that describes the evolution of a network of chemical reactions as a stochastic process. Its solution yields the probability density vector of the system at each point in time. Solving the CME numerically is in many cases computationally expensive or even infeasible as the number of reachable states can be very large or infinite. We introduce the sliding window method, which computes an approximate solution of the CME by performing a sequence of local analysis steps. In each step, only a manageable subset of states is considered, representing a "window" into the state space. In subsequent steps, the window follows the direction in which the probability mass moves, until the time period of interest has elapsed. We construct the window based on a deterministic approximation of the future behavior of the system by estimating upper and lower bounds on the populations of the chemical species.     

Efficient multipoint linkage analysis through reduction of inheritance space

Computational constraints currently limit exact multipoint linkage analysis to pedigrees of moderate size. We introduce new algorithms that allow analysis of larger pedigrees by reducing the time and memory requirements of the computation. We use the observed pedigree genotypes to reduce the number of inheritance patterns that need to be considered. The algorithms are implemented in a new version (version 2.1) of the software package GENEHUNTER. Performance gains depend on marker heterozygosity and on the number of pedigree members available for genotyping, but typically are 10-1,000-fold, compared with the performance of the previous release (version 2.0). As a result, families with up to 30 bits of inheritance information have been analyzed, and further increases in family size are feasible. In addition to computation of linkage statistics and haplotype determination, GENEHUNTER can also perform single-locus and multilocus transmission/disequilibrium tests. We describe and implement a set of permutation tests that allow determination of empirical significance levels in the presence of linkage disequilibrium among marker loci.     

Solving the chemical master equation for monomolecular reaction systems analytically

The stochastic dynamics of a well-stirred mixture of molecular species interacting through different biochemical reactions can be accurately modelled by the chemical master equation (CME). Research in the biology and scientific computing community has concentrated mostly on the development of numerical techniques to approximate the solution of the CME via many realizations of the associated Markov jump process. The domain of exact and/or efficient methods for directly solving the CME is still widely open, which is due to its large dimension that grows exponentially with the number of molecular species involved. In this article, we present an exact solution formula of the CME for arbitrary initial conditions in the case where the underlying system is governed by monomolecular reactions. The solution can be expressed in terms of the convolution of multinomial and product Poisson distributions with time-dependent parameters evolving according to the traditional reaction-rate equations. This very structured representation allows to deduce easily many properties of the solution. The model class includes many interesting examples. For more complex reaction systems, our results can be seen as a first step towards the construction of new numerical integrators, because solutions to the monomolecular case provide promising ansatz functions for Galerkin-type methods.     

Beyond the chemical master equation: Stochastic chemical kinetics coupled with auxiliary processes

The chemical master equation and its continuum approximations are indispensable tools in the modeling of chemical reaction networks. These are routinely used to capture complex nonlinear phenomena such as multimodality as well as transient events such as first-passage times, that accurately characterise a plethora of biological and chemical processes. However, some mechanisms, such as heterogeneous cellular growth or phenotypic selection at the population level, cannot be represented by the master equation and thus have been tackled separately. In this work, we propose a unifying framework that augments the chemical master equation to capture such auxiliary dynamics, and we develop and analyse a numerical solver that accurately simulates the system dynamics. We showcase these contributions by casting a diverse array of examples from the literature within this framework and applying the solver to both match and extend previous studies. Analytical calculations performed for each example validate our numerical results and benchmark the solver implementation.     

Expansion of the master equation for a biomolecular selection model

A stochastic model based on Eigen and Schuster's theory of biomolecular self-replication is studied by treating the master equation with the system-size expansion technique. The steady-state results are found to be in good agreement with the previous results and with those derived from the principle of detailed balancing. Multispecies competition and coexistence are studied carefully with the conclusions that a stable steady state is predicted for the former and a metastable state for the latter. The stochastic selection processes are also analyzed and discussed.     

DeepCME: A deep learning framework for computing solution statistics of the chemical master equation

Stochastic models of biomolecular reaction networks are commonly employed in systems and synthetic biology to study the effects of stochastic fluctuations emanating from reactions involving species with low copy-numbers. For such models, the Kolmogorov’s forward equation is called the chemical master equation (CME), and it is a fundamental system of linear ordinary differential equations (ODEs) that describes the evolution of the probability distribution of the random state-vector representing the copy-numbers of all the reacting species. The size of this system is given by the number of states that are accessible by the chemical system, and for most examples of interest this number is either very large or infinite. Moreover, approximations that reduce the size of the system by retaining only a finite number of important chemical states (e.g. those with non-negligible probability) result in high-dimensional ODE systems, even when the number of reacting species is small. Consequently, accurate numerical solution of the CME is very challenging, despite the linear nature of the underlying ODEs. One often resorts to estimating the solutions via computationally intensive stochastic simulations. The goal of the present paper is to develop a novel deep-learning approach for computing solution statistics of high-dimensional CMEs by reformulating the stochastic dynamics using Kolmogorov’s backward equation. The proposed method leverages superior approximation properties of Deep Neural Networks (DNNs) to reliably estimate expectations under the CME solution for several user-defined functions of the state-vector. This method is algorithmically based on reinforcement learning and it only requires a moderate number of stochastic simulations (in comparison to typical simulation-based approaches) to train the “policy function”. This allows not just the numerical approximation of various expectations for the CME solution but also of its sensitivities with respect to all the reaction network parameters (e.g. rate constants). We provide four examples to illustrate our methodology and provide several directions for future research.     

A master equation for a spatial population model with pair interactions

We derive a closed master equation for an individual-based population model in continuous space and time. The model and master equation include Brownian motion, reproduction via binary fission, and an interaction-dependent death rate moderated by a competition kernel. Using simulations we compare this individual-based model with the simplest approximation, the spatial logistic equation. In the limit of strong diffusion the spatial logistic equation is a good approximation to the model. However, in the limit of weak diffusion the spatial logistic equation is inaccurate because of spontaneous clustering driven by reproduction. The weak-diffusion limit can be partially analyzed using an exact solution of the master equation applicable to a competition kernel with infinite range. This analysis shows that in the case of a top-hat kernel, reducing the diffusion can increase the total population. For a Gaussian kernel, reduced diffusion invariably reduces the total population. These theoretical results are confirmed by simulation.     

A physical analysis of the ion parametric resonance model

We show, in elementary terms, using for the most part only elementary mathematics, the physical bases for the ion parametric resonance model so as to clarify the assumptions and consequences of the model. The analysis shows why, contrary to earlier conclusions, no combination of weak DC and AC magnetic fields can modify the transition rate to the ground state of excited ions. Although reinterpretations of the biological consequences of the motion of the excited ions circumvent that particular objection to the model, those changes introduce other difficulties. Also, other objections to the mechanism still stand; hence the model cannot account for any purported biological effects of weak extremely low frequency magnetic fields. Bioelectromagnetics 19:181–191, 1998. © 1998 Wiley-Liss, Inc.     

General solution of the chemical master equation and modality of marginal distributions for hierarchic first-order reaction networks

Journal of Mathematical Biology - Multimodality is a phenomenon which complicates the analysis of statistical data based exclusively on mean and variance. Here, we present criteria for...     

A genetic model of dental reduction through the probable mutation effect

A simulation approach is used in order to elucidate the nature of the hypothesized “probable mutation effect” as it applies to dental reduction in man. Computer-generated simulations of the accumulation of mutations in a human gene pool show the results of the proposed model under the influence of various parameters, as well as illustrating the nature of such genetic change through time. This approach supports a polygenic model of the probable mutation effect as a viable hypothesis for an explanation of the dental reduction which has occurred in some human populations over the last 40,000 years.     

A parametric copula model for analysis of familial binary data

Modeling the joint distribution of a binary trait (disease) within families is a tedious challenge, owing to the lack of a general statistical model with desirable properties such as the multivariate Gaussian model for a quantitative trait. Models have been proposed that either assume the existence of an underlying liability variable, the reality of which cannot be checked, or provide estimates of aggregation parameters that are dependent on the ordering of family members and on family size. We describe how a class of copula models for the analysis of exchangeable categorical data can be incorporated into a familial framework. In this class of models, the joint distribution of binary outcomes is characterized by a function of the given marginals. This function, referred to as a "copula," depends on an aggregation parameter that is weakly dependent on the marginal distributions. We propose to decompose a nuclear family into two sets of equicorrelated data (parents and offspring), each of which is characterized by an aggregation parameter (alphaFM and alphaSS, respectively). The marginal probabilities are modeled through a logistic representation. The advantage of this model is that it provides estimates of the aggregation parameters that are independent of family size and does not require any arbitrary ordering of sibs. It can be incorporated easily into segregation or combined segregation-linkage analysis and does not require extensive computer time. As an illustration, we applied this model to a combined segregation-linkage analysis of levels of plasma angiotensin I-converting enzyme (ACE) dichotomized into two classes according to the median. The conclusions of this analysis were very similar to those we had reported in an earlier familial analysis of quantitative ACE levels.     

Numerical integration of the master equation in some models of stochastic epidemiology

The processes by which disease spreads in a population of individuals are inherently stochastic. The master equation has proven to be a useful tool for modeling such processes. Unfortunately, solving the master equation analytically is possible only in limited cases (e.g., when the model is linear), and thus numerical procedures or approximation methods must be employed. Available approximation methods, such as the system size expansion method of van Kampen, may fail to provide reliable solutions, whereas current numerical approaches can induce appreciable computational cost. In this paper, we propose a new numerical technique for solving the master equation. Our method is based on a more informative stochastic process than the population process commonly used in the literature. By exploiting the structure of the master equation governing this process, we develop a novel technique for calculating the exact solution of the master equation--up to a desired precision--in certain models of stochastic epidemiology. We demonstrate the potential of our method by solving the master equation associated with the stochastic SIR epidemic model. MATLAB software that implements the methods discussed in this paper is freely available as Supporting Information S1.     

A partial differential equation model and its reduction to an ordinary differential equation model for prostate tumor growth under intermittent hormone therapy

Hormonal therapy with androgen suppression is a common treatment for advanced prostate tumors. The emergence of androgen-independent cells, however, leads to a tumor relapse under a condition of long-term androgen deprivation. Clinical trials suggest that intermittent androgen suppression (IAS) with alternating on- and off-treatment periods can delay the relapse when compared with continuous androgen suppression (CAS). In this paper, we propose a mathematical model for prostate tumor growth under IAS therapy. The model elucidates initial hormone sensitivity, an eventual relapse of a tumor under CAS therapy, and a delay of a relapse under IAS therapy, which are due to the coexistence of androgen-dependent cells, androgen-independent cells resulting from reversible changes by adaptation, and androgen-independent cells resulting from irreversible changes by genetic mutations. The model is formulated as a free boundary problem of partial differential equations that describe the evolution of populations of the abovementioned three types of cells during on-treatment periods and off-treatment periods. Moreover, the model can be transformed into a piecewise linear ordinary differential equation model by introducing three new volume variables, and the study of the resulting model may help to devise optimal IAS schedules.     

Reorganising the order Bacillales through phylogenomics

Bacterial classification at higher taxonomic ranks such as the order and family levels is currently reliant on phylogenetic analysis of 16S rRNA and the presence of shared phenotypic characteristics. However, these may not be reflective of the true genotypic and phenotypic relationships of taxa. This is evident in the order Bacillales, members of which are defined as aerobic, spore-forming and rod-shaped bacteria. However, some taxa are anaerobic, asporogenic and coccoid. 16S rRNA gene phylogeny is also unable to elucidate the taxonomic positions of several families incertae sedis within this order. Whole genome-based phylogenetic approaches may provide a more accurate means to resolve higher taxonomic levels. A suite of phylogenomic approaches were applied to re-evaluate the taxonomy of 80 representative taxa of eight families (and six family incertae sedis taxa) within the order Bacillales. This showed several anomalies in the current family and order level classifications including the existence of four Bacillaceae and two Paenibacillaceae “family” clades. Furthermore, the families Staphylococcaceae and Listeriaceae belong to the sister order Lactobacillales. Finally, we propose a consensus phylogenomic approach which may diminish algorithmic biases associated with single approaches and facilitate more accurate classification of a broad range of taxa at the higher taxonomic levels.