Analysis of stochastic strategies in bacterial competence: a master equation approach

Competence is a transiently differentiated state that certain bacterial cells reach when faced with a stressful environment. Entrance into competence can be attributed to the excitability of the dynamics governing the genetic circuit that regulates this cellular behavior. Like many biological behaviors, entrance into competence is a stochastic event. In this case cellular noise is responsible for driving the cell from a vegetative state into competence and back. In this work we present a novel numerical method for the analysis of stochastic biochemical events and use it to study the excitable dynamics responsible for competence in Bacillus subtilis. Starting with a Finite State Projection (FSP) solution of the chemical master equation (CME), we develop efficient numerical tools for accurately computing competence probability. Additionally, we propose a new approach for the sensitivity analysis of stochastic events and utilize it to elucidate the robustness properties of the competence regulatory genetic circuit. We also propose and implement a numerical method to calculate the expected time it takes a cell to return from competence. Although this study is focused on an example of cell-differentiation in Bacillus subtilis, our approach can be applied to a wide range of stochastic phenomena in biological systems.     

Extracellular noise-induced stochastic synchronization in heterogeneous quorum sensing network

Quorum sensing is a bacterial mechanism used to synchronize the coordinated response of a microbial population. Because quorum sensing in Gram-negative bacteria depends on release and detection of a diffusible signaling molecule (autoinducer) among a multicellular group, it is considered a simple form of cell-cell communication for the purposes of mathematical analysis. Stochastic equation systems have provided a common approach to model biochemical or biophysical processes. Recently, the effect of noise to synchronize a specific homogeneous quorum sensing network was successfully modeled using a stochastic equation system with fixed parameters. The question remains of how to model quorum sensing networks in a general setting. To address this question, we first set a stochastic equation system as a general model for a heterogeneous quorum sensing network. Then, using two relevant biophysical characteristics of Gram-negative bacteria (the permeability of the cell membrane to the autoinducer and the symmetry of autoinducer diffusion) we construct the solution of the stochastic equation system at an abstract level. The solution indicates that stable synchronization of a quorum sensing network is robustly induced by an environment with a heterogenous distribution of extracellular and intracellular noise. The synchronization is independent of the initial state of the system and is solely the result of the connectivity of the cell network established through the effects of extracellular noise.     

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.     

Investigating the two-moment characterisation of subcellular biochemical networks

While ordinary differential equations (ODEs) form the conceptual framework for modelling many cellular processes, specific situations demand stochastic models to capture the influence of noise. The most common formulation of stochastic models for biochemical networks is the chemical master equation (CME). While stochastic simulations are a practical way to realise the CME, analytical approximations offer more insight into the influence of noise. Towards that end, the two-moment approximation (2MA) is a promising addition to the established analytical approaches including the chemical Langevin equation (CLE) and the related linear noise approximation (LNA). The 2MA approach directly tracks the mean and (co)variance which are coupled in general. This coupling is not obvious in CME and CLE and ignored by LNA and conventional ODE models. We extend previous derivations of 2MA by allowing (a) non-elementary reactions and (b) relative concentrations. Often, several elementary reactions are approximated by a single step. Furthermore, practical situations often require the use of relative concentrations. We investigate the applicability of the 2MA approach to the well-established fission yeast cell cycle model. Our analytical model reproduces the clustering of cycle times observed in experiments. This is explained through multiple resettings of M-phase promoting factor (MPF), caused by the coupling between mean and (co)variance, near the G2/M transition.     

Current noise around steady states in discrete transport systems

Subject of this paper is the transport noise in discrete systems. The transport systems are given by a number (n) of binding sites separated by energy barriers. These binding sites may be in contact with constant outer reservoirs. The state of the system is characterized by the occupation numbers of particles (current carriers) at these binding sites. The change in time of the occupation numbers is generated by individual “jumps” of particles over the energy barriers, building up the flux matter (for charged particles: the electric current). In the limit n → ∞ continuum processes as e.g. usual diffusion are included in the transport model. The fluctuations in occupation numbers and other quantities linearly coupled to the occupation numbers may be treated with the usual master equation approach. The treatment of the fluctuations in fluxes (current) makes necessary a different theoretical approach which is presented in this paper under the assumption of vanishing interactions between the particles. This approach may be applied to a number of different transport systems in biology and physics (ion transport through porous channels in membranes, carrier mediated ion transport through membranes, jump diffusion e.g. in superionic conductors). As in the master equation approach the calculation of correlations and noise spectra may be reduced to the solution of the macroscopic equations for the occupation numbers. This result may be regarded as a generalization to non-equilibrium current fluctuations of the usual Nyquist theorem relating the current (voltage) noise spectrum in thermal equilibrium to the macroscopic frequency dependent admittance.The validity of the general approach is demonstrated by the calculation of the autocorrelation function and spectrum of current noise for a number of special examples (e.g, pores in membrances, carrier mediated ion transport).     

Noise suppression in stochastic genetic circuits using PID controllers

Inside individual cells, protein population counts are subject to molecular noise due to low copy numbers and the inherent probabilistic nature of biochemical processes. We investigate the effectiveness of proportional, integral and derivative (PID) based feedback controllers to suppress protein count fluctuations originating from two noise sources: bursty expression of the protein, and external disturbance in protein synthesis. Designs of biochemical reactions that function as PID controllers are discussed, with particular focus on individual controllers separately, and the corresponding closed-loop system is analyzed for stochastic controller realizations. Our results show that proportional controllers are effective in buffering protein copy number fluctuations from both noise sources, but this noise suppression comes at the cost of reduced static sensitivity of the output to the input signal. In contrast, integral feedback has no effect on the protein noise level from stochastic expression, but significantly minimizes the impact of external disturbances, particularly when the disturbance comes at low frequencies. Counter-intuitively, integral feedback is found to amplify external disturbances at intermediate frequencies. Next, we discuss the design of a coupled feedforward-feedback biochemical circuit that approximately functions as a derivate controller. Analysis using both analytical methods and Monte Carlo simulations reveals that this derivative controller effectively buffers output fluctuations from bursty stochastic expression, while maintaining the static input-output sensitivity of the open-loop system. In summary, this study provides a systematic stochastic analysis of biochemical controllers, and paves the way for their synthetic design and implementation to minimize deleterious fluctuations in gene product levels.     

Formulas for intrinsic noise evaluation in oscillatory genetic networks

The linear noise approximation is a useful method for stochastic noise evaluations in genetic regulatory networks, where the covariance equation described as a Lyapunov equation plays a central role. We discuss the linear noise approximation method for evaluations of an intrinsic noise in autonomously oscillatory genetic networks; in such oscillatory networks, the covariance equation becomes a periodic differential equation that provides generally an unbounded covariance matrix, so that the standard method of noise evaluation based on the covariance matrix cannot be adopted directly. In this paper, we develop a new method of noise evaluation in oscillatory genetic networks; first, we investigate structural properties, e.g., orbital stability and periodicity, of the solutions to the covariance equation given as a periodic Lyapunov differential equation by using the Floquet-Lyapunov theory, and propose a global measure for evaluating stochastic amplitude fluctuations on the periodic trajectory; we also derive an evaluation formula for the period fluctuation. Finally, we apply our method to a model of circadian oscillations based on negative auto-regulation of gene expression, and show validity of our method by comparing the evaluation results with stochastic simulations.     

Noise-induced modulation of the relaxation kinetics around a non-equilibrium steady state of non-linear chemical reaction networks

Stochastic effects from correlated noise non-trivially modulate the kinetics of non-linear chemical reaction networks. This is especially important in systems where reactions are confined to small volumes and reactants are delivered in bursts. We characterise how the two noise sources confinement and burst modulate the relaxation kinetics of a non-linear reaction network around a non-equilibrium steady state. We find that the lifetimes of species change with burst input and confinement. Confinement increases the lifetimes of all species that are involved in any non-linear reaction as a reactant. Burst monotonically increases or decreases lifetimes. Competition between burst-induced and confinement-induced modulation may hence lead to a non-monotonic modulation. We quantify lifetime as the integral of the time autocorrelation function (ACF) of concentration fluctuations around a non-equilibrium steady state of the reaction network. Furthermore, we look at the first and second derivatives of the ACF, each of which is affected in opposite ways by burst and confinement. This allows discriminating between these two noise sources. We analytically derive the ACF from the linear Fokker-Planck approximation of the chemical master equation in order to establish a baseline for the burst-induced modulation at low confinement. Effects of higher confinement are then studied using a partial-propensity stochastic simulation algorithm. The results presented here may help understand the mechanisms that deviate stochastic kinetics from its deterministic counterpart. In addition, they may be instrumental when using fluorescence-lifetime imaging microscopy (FLIM) or fluorescence-correlation spectroscopy (FCS) to measure confinement and burst in systems with known reaction rates, or, alternatively, to correct for the effects of confinement and burst when experimentally measuring reaction rates.     

Stochastic Oscillations Induced by Intrinsic Fluctuations in a Self-Repressing Gene

Biochemical reaction networks are subjected to large fluctuations attributable to small molecule numbers, yet underlie reliable biological functions. Thus, it is important to understand how regularity can emerge from noise. Here, we study the stochastic dynamics of a self-repressing gene with arbitrarily long or short response time. We find that when the mRNA and protein half-lives are approximately equal to the gene response time, fluctuations can induce relatively regular oscillations in the protein concentration. To gain insight into this phenomenon at the crossroads of determinism and stochasticity, we use an intermediate theoretical approach, based on a moment-closure approximation of the master equation, which allows us to take into account the binary character of gene activity. We thereby obtain differential equations that describe how nonlinearity can feed-back fluctuations into the mean-field equations to trigger oscillations. Finally, our results suggest that the self-repressing Hes1 gene circuit exploits this phenomenon to generate robust oscillations, inasmuch as its time constants satisfy precisely the conditions we have identified.     

Cellular growth and division in the Gillespie algorithm

Recent experimental studies elucidating the importance of noise in gene regulation have ignited widespread interest in Gillespie's stochastic simulation technique for biochemical networks. We formulate modifications to the Gillespie algorithm which are necessary to correctly simulate chemical reactions with time-dependent reaction rates. We concentrate on time dependence of kinetic rates arising from the periodic process of growth and division of the cellular volume, and demonstrate that a careful re-derivation of the Gillespie algorithm is important when all stochastically simulated reactions have rates slower or comparable to the cellular growth rate. For an unregulated single-gene system, we illustrate our findings using recently proposed hybrid simulation techniques, and systematically compare our algorithm with analytic results obtained from the chemical master equation.     

Moment fitting for parameter inference in repeatedly and partially observed stochastic biological models

The inference of reaction rate parameters in biochemical network models from time series concentration data is a central task in computational systems biology. Under the assumption of well mixed conditions the network dynamics are typically described by the chemical master equation, the Fokker Planck equation, the linear noise approximation or the macroscopic rate equation. The inverse problem of estimating the parameters of the underlying network model can be approached in deterministic and stochastic ways, and available methods often compare individual or mean concentration traces obtained from experiments with theoretical model predictions when maximizing likelihoods, minimizing regularized least squares functionals, approximating posterior distributions or sequentially processing the data. In this article we assume that the biological reaction network can be observed at least partially and repeatedly over time such that sample moments of species molecule numbers for various time points can be calculated from the data. Based on the chemical master equation we furthermore derive closed systems of parameter dependent nonlinear ordinary differential equations that predict the time evolution of the statistical moments. For inferring the reaction rate parameters we suggest to not only compare the sample mean with the theoretical mean prediction but also to take the residual of higher order moments explicitly into account. Cost functions that involve residuals of higher order moments may form landscapes in the parameter space that have more pronounced curvatures at the minimizer and hence may weaken or even overcome parameter sloppiness and uncertainty. As a consequence both deterministic and stochastic parameter inference algorithms may be improved with respect to accuracy and efficiency. We demonstrate the potential of moment fitting for parameter inference by means of illustrative stochastic biological models from the literature and address topics for future research.     

Stochastic modeling of nonlinear epidemiology

The objectives of this paper to analyse, model and simulate the spread of an infectious disease by resorting to modern stochastic algorithms. The approach renders it possible to circumvent the simplifying assumption of linearity imposed in the majority of the past works on stochastic analysis of epidemic processes. Infectious diseases are often transmitted through contacts of those infected with those susceptible; hence the processes are inherently nonlinear. According to the classical model of Kermack and McKendrick, or the SIR model, three classes of populations are involved in two types of processes: conversion of susceptibles (S) to infectives (I) and conversion of infectives to removed (R). The master equations of the SIR process have been formulated through the probabilistic population balance around a particular state by considering the mutually exclusive events. The efficacy of the present methodology is mainly attributable to its ability to derive the governing equations for the means, variances and covariance of the random variables by the method of system-size expansion of the nonlinear master equations. Solving these equations simultaneously along with rates associated influenza epidemic data yields information concerning not only the means of the three populations but also the minimal uncertainties of these populations inherent in the epidemic. The stochastic pathways of the three different classes of populations during an epidemic, i.e. their means and the fluctuations around these means, have also been numerically simulated independently by the algorithm derived from the master equations, as well as by an event-driven Monte Carlo algorithm. The master equation and Monte Carlo algorithms have given rise to the identical results.     

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 decision-making Fokker–Planck model in computational neuroscience

In computational neuroscience, decision-making may be explained analyzing models based on the evolution of the average firing rates of two interacting neuron populations, e.g., in bistable visual perception problems. These models typically lead to a multi-stable scenario for the concerned dynamical systems. Nevertheless, noise is an important feature of the model taking into account both the finite-size effects and the decision’s robustness. These stochastic dynamical systems can be analyzed by studying carefully their associated Fokker–Planck partial differential equation. In particular, in the Fokker–Planck setting, we analytically discuss the asymptotic behavior for large times towards a unique probability distribution, and we propose a numerical scheme capturing this convergence. These simulations are used to validate deterministic moment methods recently applied to the stochastic differential system. Further, proving the existence, positivity and uniqueness of the probability density solution for the stationary equation, as well as for the time evolving problem, we show that this stabilization does happen. Finally, we discuss the convergence of the solution for large times to the stationary state. Our approach leads to a more detailed analytical and numerical study of decision-making models applied in computational neuroscience.     

Non-equilibrium voltage noise generated by ion transport through pores

In this paper, we describe a systematic approach to the theoretical analysis of non-equilibrium voltage noise that arises from ions moving through pores in membranes. We assume that an ion must cross one or two barriers in the pore in order to move from one side of the membrane to the other. In our analysis, we consider the following factors: a) surface charge as a variable in the kinetic equations, b) linearization of the kinetic equations, c) master equation approach to fluctuations. To analyze the voltage noise arising from ion movement through a two barrier (i.e., one binding site) pore, we included the effects of ions in the channel's interior on the voltage noise. The current clamp is considered as a white noise generating additional noise in the system. In contrast to what is found for current noise, at low frequencies the voltage noise intensity is reduced by increasing voltage across the membrane. With this approach, we demonstrate explicity for the examples treated that, apart from additional noise generated by the current clamp, the non-equilibrium voltage fluctuations can be related to the current fluctuations by the complex admittance.     

Enzymes as molecular automata: a stochastic model of self-oscillatory glycolytic cycles in cellular metabolism

A stochastic model based on the molecular automata approach was developed to simulate the cyclic dynamics of glycolysis-gluconeogenesis in cell energy metabolism. The stochastic algorithm, based on Gillespie's method, simulates the master equation associated with any network of enzymatically controlled reactions. This model of the glycolytic-gluconeogenetic cycle was developed by assembling the stochastic kinetic reactions controlled by two enzymes: phosphofructokinase (PFKase) and fructose-1, 6-biphosphatase (FBPase). In order to obtain the hysteresis behaviour predicted by classical Sel'kov analysis, a minimum complexity is required in the allosteric regulations. This implies that the FBPase cannot have a single binding site for its transition to the inactive state (a minimum of two or three binding sites is necessary). Given the multimeric structure of this enzyme, this kinetic hypothesis is tenable. Other possible applications of the stochastic automata approach for different cases of channels, receptors and enzymatic control are suggested.