Extended Kalman Filter for Estimation of Parameters in Nonlinear State-Space Models of Biochemical Networks

It is system dynamics that determines the function of cells, tissues and organisms. To develop mathematical models and estimate their parameters are an essential issue for studying dynamic behaviors of biological systems which include metabolic networks, genetic regulatory networks and signal transduction pathways, under perturbation of external stimuli. In general, biological dynamic systems are partially observed. Therefore, a natural way to model dynamic biological systems is to employ nonlinear state-space equations. Although statistical methods for parameter estimation of linear models in biological dynamic systems have been developed intensively in the recent years, the estimation of both states and parameters of nonlinear dynamic systems remains a challenging task. In this report, we apply extended Kalman Filter (EKF) to the estimation of both states and parameters of nonlinear state-space models. To evaluate the performance of the EKF for parameter estimation, we apply the EKF to a simulation dataset and two real datasets: JAK-STAT signal transduction pathway and Ras/Raf/MEK/ERK signaling transduction pathways datasets. The preliminary results show that EKF can accurately estimate the parameters and predict states in nonlinear state-space equations for modeling dynamic biochemical networks.     

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.     

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.     

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.     

Direct Solution of the Chemical Master Equation Using Quantized Tensor Trains

The Chemical Master Equation (CME) is a cornerstone of stochastic analysis and simulation of models of biochemical reaction networks. Yet direct solutions of the CME have remained elusive. Although several approaches overcome the infinite dimensional nature of the CME through projections or other means, a common feature of proposed approaches is their susceptibility to the curse of dimensionality, i.e. the exponential growth in memory and computational requirements in the number of problem dimensions. We present a novel approach that has the potential to “lift” this curse of dimensionality. The approach is based on the use of the recently proposed Quantized Tensor Train (QTT) formatted numerical linear algebra for the low parametric, numerical representation of tensors. The QTT decomposition admits both, algorithms for basic tensor arithmetics with complexity scaling linearly in the dimension (number of species) and sub-linearly in the mode size (maximum copy number), and a numerical tensor rounding procedure which is stable and quasi-optimal. We show how the CME can be represented in QTT format, then use the exponentially-converging -discontinuous Galerkin discretization in time to reduce the CME evolution problem to a set of QTT-structured linear equations to be solved at each time step using an algorithm based on Density Matrix Renormalization Group (DMRG) methods from quantum chemistry. Our method automatically adapts the “basis” of the solution at every time step guaranteeing that it is large enough to capture the dynamics of interest but no larger than necessary, as this would increase the computational complexity. Our approach is demonstrated by applying it to three different examples from systems biology: independent birth-death process, an example of enzymatic futile cycle, and a stochastic switch model. The numerical results on these examples demonstrate that the proposed QTT method achieves dramatic speedups and several orders of magnitude storage savings over direct approaches.     

Identification of dynamic mass-action biochemical reaction networks using sparse Bayesian methods

Identifying the reactions that govern a dynamical biological system is a crucial but challenging task in systems biology. In this work, we present a data-driven method to infer the underlying biochemical reaction system governing a set of observed species concentrations over time. We formulate the problem as a regression over a large, but limited, mass-action constrained reaction space and utilize sparse Bayesian inference via the regularized horseshoe prior to produce robust, interpretable biochemical reaction networks, along with uncertainty estimates of parameters. The resulting systems of chemical reactions and posteriors inform the biologist of potentially several reaction systems that can be further investigated. We demonstrate the method on two examples of recovering the dynamics of an unknown reaction system, to illustrate the benefits of improved accuracy and information obtained.     

Analysis of operating principles with S-system models

Operating principles address general questions regarding the response dynamics of biological systems as we observe or hypothesize them, in comparison to a priori equally valid alternatives. In analogy to design principles, the question arises: Why are some operating strategies encountered more frequently than others and in what sense might they be superior? It is at this point impossible to study operation principles in complete generality, but the work here discusses the important situation where a biological system must shift operation from its normal steady state to a new steady state. This situation is quite common and includes many stress responses. We present two distinct methods for determining different solutions to this task of achieving a new target steady state. Both methods utilize the property of S-system models within Biochemical Systems Theory (BST) that steady states can be explicitly represented as systems of linear algebraic equations. The first method uses matrix inversion, a pseudo-inverse, or regression to characterize the entire admissible solution space. Operations on the basis of the solution space permit modest alterations of the transients toward the target steady state. The second method uses standard or mixed integer linear programming to determine admissible solutions that satisfy criteria of functional effectiveness, which are specified beforehand. As an illustration, we use both methods to characterize alternative response patterns of yeast subjected to heat stress, and compare them with observations from the literature.     

Parameter estimation in stochastic biochemical reactions

Gene regulatory, signal transduction and metabolic networks are major areas of interest in the newly emerging field of systems biology. In living cells, stochastic dynamics play an important role; however, the kinetic parameters of biochemical reactions necessary for modelling these processes are often not accessible directly through experiments. The problem of estimating stochastic reaction constants from molecule count data measured, with error, at discrete time points is considered. For modelling the system, a hidden Markov process is used, where the hidden states are the true molecule counts, and the transitions between those states correspond to reaction events following collisions of molecules. Two different algorithms are proposed for estimating the unknown model parameters. The first is an approximate maximum likelihood method that gives good estimates of the reaction parameters in systems with few possible reactions in each sampling interval. The second algorithm, treating the data as exact measurements, approximates the number of reactions in each sampling interval by solving a simple linear equation. Maximising the likelihood based on these approximations can provide good results, even in complex reaction systems.     

Time-dependent propagators for stochastic models of gene expression: an analytical method

The inherent stochasticity of gene expression in the context of regulatory networks profoundly influences the dynamics of the involved species. Mathematically speaking, the propagators which describe the evolution of such networks in time are typically defined as solutions of the corresponding chemical master equation (CME). However, it is not possible in general to obtain exact solutions to the CME in closed form, which is due largely to its high dimensionality. In the present article, we propose an analytical method for the efficient approximation of these propagators. We illustrate our method on the basis of two categories of stochastic models for gene expression that have been discussed in the literature. The requisite procedure consists of three steps: a probability-generating function is introduced which transforms the CME into (a system of) partial differential equations (PDEs); application of the method of characteristics then yields (a system of) ordinary differential equations (ODEs) which can be solved using dynamical systems techniques, giving closed-form expressions for the generating function; finally, propagator probabilities can be reconstructed numerically from these expressions via the Cauchy integral formula. The resulting ‘library’ of propagators lends itself naturally to implementation in a Bayesian parameter inference scheme, and can be generalised systematically to related categories of stochastic models beyond the ones considered here.

     

Variable elimination in post-translational modification reaction networks with mass-action kinetics

We define a subclass of chemical reaction networks called post-translational modification systems. Important biological examples of such systems include MAPK cascades and two-component systems which are well-studied experimentally as well as theoretically. The steady states of such a system are solutions to a system of polynomial equations. Even for small systems the task of finding the solutions is daunting. We develop a mathematical framework based on the notion of a cut (a particular subset of species in the system), which provides a linear elimination procedure to reduce the number of variables in the system to a set of core variables. The steady states are parameterized algebraically by the core variables, and graphical conditions for when steady states with positive core variables imply positivity of all variables are given. Further, minimal cuts are the connected components of the species graph and provide conservation laws. A criterion for when a (maximal) set of independent conservation laws can be derived from cuts is given.     

The numerical solution of stiff differential equations

This paper first discusses the conditions in which a set of differential equations should give stable solutions, starting with linear systems assuming that these do not differ greatly in this respect from non-linear systems. Methods of investigating the stability of particular systems are briefly discussed. Most real biochemical systems are known from observation to be stable, but little is known of the regions over which stability persists; moreover, models of biochemical systems may not be stable, because of inaccurate choice of parameter values.The separate problem of stability and accuracy in numerical methods of approximating the solution of systems of non-linear equations is then treated. Stress is laid on the consistently unsatisfactory results given by explicit methods for systems containing "stiff" equations, and implicit multistep methods are particularly recommended for this class of problem, which is likely to include many biochemical model systems. Finally, an iteration procedure likely to give convergence both in multistep methods and in the steady-state approach is recommended, and areas in which improvement in methods is likely to occur are outlined.     

Explaining Differences in Repatriation Experiences: The Discovery of Coupled and Decoupled Systems

We report the results of a four-year, multiphase study on the overseas assignment at General Motors Corporation (GM). Our objective is to explain variation in the repatriation experiences of International Service Personnel (ISPs). While our principal focus is the documentation of an inductive discovery process facilitated through interviews with a wide variety of GM employees, we also present a preliminary test of an explanatory hypothesis. The hypothesis states that variation in repatriation experiences is related directly to variation in the structure and ideology of organizational units that send and receive ISPs. Units with direct linkages between GM's domestic and overseas operations (coupled systems) are associated with a pro-international ideology and positive ISP repatriation experiences. Units with no administrative or operational linkages between the domestic and overseas arenas (decoupled systems) are associated with an anti-international ideology and negative ISP repatriation experiences. Our findings suggest that the structural and ideological properties of organizational units affect both repatriation from overseas assignments and ISP career paths.     

PathMiner: predicting metabolic pathways by heuristic search

MOTIVATION: Automated methods for biochemical pathway inference are becoming increasingly important for understanding biological processes in living and synthetic systems. With the availability of data on complete genomes and increasing information about enzyme-catalyzed biochemistry it is becoming feasible to approach this problem computationally. In this paper we present PathMiner, a system for automatic metabolic pathway inference. PathMiner predicts metabolic routes by reasoning over transformations using chemical and biological information. RESULTS: We build a biochemical state-space using data from known enzyme-catalyzed transformations in Ligand, including, 2917 unique transformations between 3890 different compounds. To predict metabolic pathways we explore this state-space by developing an informed search algorithm. For this purpose we develop a chemically motivated heuristic to guide the search. Since the algorithm does not depend on predefined pathways, it can efficiently identify plausible routes using known biochemical transformations.     

Aging-related changes of EEG synchronization during a visual working memory task

Differences of EEG synchronization between normal old and young people during a working memory (WM) task were investigated. The synchronization likelihood (SL) is a novel method to assessed synchronization in multivariate time series for non-stationary systems. To evaluate this method to study the mechanisms of WM, we calculated the SL values in brain electrical activity for both resting state and task state. EEG signals were recorded from 14 young adults and 12 old adults during two different states, respectively. SL was used to measure EEG synchronization between 19 electrodes in delta, theta, alpha1, alpha2 and beta frequency bands. Bad task performance and significantly decreased EEG synchronization were found in old group compared to young group in alpha1, alpha2 and beta frequency bands during the WM task. Moreover, significantly decreased EEG synchronization in beta band in the elder was also detected during the resting state. The findings suggested that reduced EEG synchronization may be one of causes for WM capacity decline along with healthy aging.     

On the stationary state analysis of reaction-diffusion mechanisms for biological pattern formation

We present a biologically plausible two-variable reaction-diffusion model for the developing vertebrate limb, for which we postulate the existence of a stationary solution. A consequence of this assumption is that the stationary state depends on only a single concentration-variable. Under these circumstances, features of potential biological significance, such as the dependence of the steady-state concentration profile of this variable on parameters such as tissue size and shape, can be studied without detailed information about the rate functions. As the existence and stability of stationary solutions, which must be assumed for any biochemical system governing morphogenesis, cannot be investigated without such information, an analysis is made of the minimal requirements for stable, stationary non-uniform solutions in a general class of reaction-diffusion systems. We discuss the strategy of studying stationary-state properties of systems that are incompletely specified. Where abrupt transitions between successive compartment-sizes occur, as in the developing limb, we argue that it is reasonable to model pattern reorganization as a sequence of independent stationary states.     

Determination of the quaternary phase diagram of the water-ethylene glycol-sucrose-NaCl system and a comparison between two theoretical methods for synthetic phase diagrams

Characterization of the thermodynamic properties of multi-solute aqueous solutions is of critical importance for biological and biochemical research. For example, the phase diagrams of aqueous systems, containing salts, saccharides, and plasma membrane permeating solutes, are indispensible in the field of cryobiology and pharmacology. However, only a few ternary phase diagrams are currently available for these systems. In this study, an auto-sampler differential scanning calorimeter (DSC) was used to determine the quaternary phase diagram of the water-ethylene glycol-sucrose-NaCl system. To improve the accuracy of melting point measurement, a “mass-redemption” method was also applied for the DSC technique. Base on the analyses of these experimental data, a comparison was made between the two practical approaches to generate phase diagrams of multi-solute solutions from those of single-solute solutions: the summation of cubic polynomial melting point equations versus the use of osmotic virial equations with cross coefficients. The calculated values of the model standard deviations suggested that both methods are satisfactory for characterizing this quaternary system.     

A Hybrid of the Chemical Master Equation and the Gillespie Algorithm for Efficient Stochastic Simulations of Sub-Networks

Modeling stochastic behavior of chemical reaction networks is an important endeavor in many aspects of chemistry and systems biology. The chemical master equation (CME) and the Gillespie algorithm (GA) are the two most fundamental approaches to such modeling; however, each of them has its own limitations: the GA may require long computing times, while the CME may demand unrealistic memory storage capacity. We propose a method that combines the CME and the GA that allows one to simulate stochastically a part of a reaction network. First, a reaction network is divided into two parts. The first part is simulated via the GA, while the solution of the CME for the second part is fed into the GA in order to update its propensities. The advantage of this method is that it avoids the need to solve the CME or stochastically simulate the entire network, which makes it highly efficient. One of its drawbacks, however, is that most of the information about the second part of the network is lost in the process. Therefore, this method is most useful when only partial information about a reaction network is needed. We tested this method against the GA on two systems of interest in biology - the gene switch and the Griffith model of a genetic oscillator—and have shown it to be highly accurate. Comparing this method to four different stochastic algorithms revealed it to be at least an order of magnitude faster than the fastest among them.     

Inference of biological S-system using the separable estimation method and the genetic algorithm

Reconstruction of a biological system from its experimental time series data is a challenging task in systems biology. The S-system which consists of a group of nonlinear ordinary differential equations (ODEs) is an effective model to characterize molecular biological systems and analyze the system dynamics. However, inference of S-systems without the knowledge of system structure is not a trivial task due to its nonlinearity and complexity. In this paper, a pruning separable parameter estimation algorithm (PSPEA) is proposed for inferring S-systems. This novel algorithm combines the separable parameter estimation method (SPEM) and a pruning strategy, which includes adding an l? regularization term to the objective function and pruning the solution with a threshold value. Then, this algorithm is combined with the continuous genetic algorithm (CGA) to form a hybrid algorithm that owns the properties of these two combined algorithms. The performance of the pruning strategy in the proposed algorithm is evaluated from two aspects: the parameter estimation error and structure identification accuracy. The results show that the proposed algorithm with the pruning strategy has much lower estimation error and much higher identification accuracy than the existing method.     

Protein import into yeast mitochondria: the inner membrane import site protein ISP45 is the MPI1 gene product.

Protein import across both mitochondrial membranes is mediated by the cooperation of two distinct protein transport systems, one in the outer and the other in the inner membrane. Previously we described a 45 kDa yeast mitochondrial inner membrane protein (ISP45) that can be cross-linked to a partially translocated precursor protein (Scherer et al., 1992). We have now purified ISP45 to homogeneity and identified it as the product of the nuclear MPI1 gene. Identity of ISP45 with the MPI1 gene product was shown by microsequencing of three tryptic ISP45 peptides and by demonstrating that an antibody against an Mpi1p-beta-galactosidase fusion protein specifically recognizes ISP45. Antibodies monospecific for ISP45 inhibited protein import into right-side-out mitochondrial inner membrane vesicles, but not into intact mitochondria. On solubilizing mitochondria, ISP45 was rapidly converted to a 40 kDa proteolytic fragment unless mitochondria were first denatured with trichloroacetic acid. The combined genetic and biochemical evidence identifies ISP45/Mpi1p as a component of the protein import system of the yeast mitochondrial inner membrane.