Approximate method for stochastic chemical kinetics with two-time scales by chemical Langevin equations.

The frequently used reduction technique is based on the chemical master equation for stochastic chemical kinetics with two-time scales, which yields the modified stochastic simulation algorithm (SSA). For the chemical reaction processes involving a large number of molecular species and reactions, the collection of slow reactions may still include a large number of molecular species and reactions. Consequently, the SSA is still computationally expensive. Because the chemical Langevin equations (CLEs) can effectively work for a large number of molecular species and reactions, this paper develops a reduction method based on the CLE by the stochastic averaging principle developed in the work of Khasminskii and Yin [SIAM J. Appl. Math. 56, 1766-1793 (1996); ibid. 56, 1794-1819 (1996)] to average out the fast-reacting variables. This reduction method leads to a limit averaging system, which is an approximation of the slow reactions. Because in the stochastic chemical kinetics, the CLE is seen as the approximation of the SSA, the limit averaging system can be treated as the approximation of the slow reactions. As an application, we examine the reduction of computation complexity for the gene regulatory networks with two-time scales driven by intrinsic noise. For linear and nonlinear protein production functions, the simulations show that the sample average (expectation) of the limit averaging system is close to that of the slow-reaction process based on the SSA. It demonstrates that the limit averaging system is an efficient approximation of the slow-reaction process in the sense of the weak convergence.

Parameter estimation in stochastic chemical kinetic models using derivative free optimization and bootstrapping.

Recent years have seen increasing popularity of stochastic chemical kinetic models due to their ability to explain and model several critical biological phenomena. Several developments in high resolution fluorescence microscopy have enabled researchers to obtain protein and mRNA data on the single cell level. The availability of these data along with the knowledge that the system is governed by a stochastic chemical kinetic model leads to the problem of parameter estimation. This paper develops a new method of parameter estimation for stochastic chemical kinetic models. There are three components of the new method. First, we propose a new expression for likelihood of the experimental data. Second, we use sample path optimization along with UOBYQA-Fit, a variant of of Powell's unconstrained optimization by quadratic approximation, for optimization. Third, we use a variant of Efron's percentile bootstrapping method to estimate the confidence regions for the parameter estimates. We apply the parameter estimation method in an RNA dynamics model of E. coli. We test the parameter estimates obtained and the confidence regions in this model. The testing of the parameter estimation method demonstrates the efficiency, reliability, and accuracy of the new method.

Comparison of finite difference based methods to obtain sensitivities of stochastic chemical kinetic models.

Sensitivity analysis is a powerful tool in determining parameters to which the system output is most responsive, in assessing robustness of the system to extreme circumstances or unusual environmental conditions, in identifying rate limiting pathways as a candidate for drug delivery, and in parameter estimation for calculating the Hessian of the objective function. Anderson [SIAM J. Numer. Anal. 50, 2237 (2012)] shows the advantages of the newly developed coupled finite difference (CFD) estimator over the common reaction path (CRP) [M. Rathinam, P. W. Sheppard, and M. Khammash, J. Chem. Phys. 132, 034103 (2010)] estimator. In this paper, we demonstrate the superiority of the CFD estimator over the common random number (CRN) estimator in a number of scenarios not considered previously in the literature, including the sensitivity of a negative log likelihood function for parameter estimation, the sensitivity of being in a rare state, and a sensitivity with fast fluctuating species. In all examples considered, the superiority of CFD over CRN is demonstrated. We also provide an example in which the CRN method is superior to the CRP method, something not previously observed in the literature. These examples, along with Anderson's results, lead to the conclusion that CFD is currently the best estimator in the class of finite difference estimators of stochastic chemical kinetic models.

Markedly Improved Catalytic Dehydration of Sorbitol to Isosorbide by Sol-Gel Sulfated Zirconia: A Quantitative Structure-Reactivity Study.

Isosorbide, a bicyclic C6 diol, has considerable value as a precursor for the production of bio-derived polymers. Current production of isosorbide from sorbitol utilizes homogeneous acid, commonly H2SO4, creating harmful waste and complicating separation. Thus, a heterogeneous acid catalyst capable of producing isosorbide from sorbitol in high yield under mild conditions would be desirable. Reported here is a quantitative investigation of the liquid-phase dehydration of neat sorbitol over sulfated zirconia (SZ) solid acid catalysts produced via sol-gel synthesis. The catalyst preparation allows for precise surface area control (morphology) and tunable catalytic activity. The S/Zr ratio (0.1-2.0) and calcination temperature (425-625 °C) were varied to evaluate their effects on morphology, acidity, and reaction kinetics and, as a result, to optimize the catalytic system for this transformation. With the optimal SZ catalyst, complete conversion of sorbitol occurred in <2 h under mild conditions to give isosorbide in 76% yield. Overall, the quantitative kinetics and structure-reactivity studies provided valuable insights into the parameters that govern product yields and SZ catalyst activity, central among these being the relative proportion of acid site types and Brønsted surface density.

The stochastic quasi-steady-state assumption: reducing the model but not the noise.

Highly reactive species at small copy numbers play an important role in many biological reaction networks. We have described previously how these species can be removed from reaction networks using stochastic quasi-steady-state singular perturbation analysis (sQSPA). In this paper we apply sQSPA to three published biological models: the pap operon regulation, a biochemical oscillator, and an intracellular viral infection. These examples demonstrate three different potential benefits of sQSPA. First, rare state probabilities can be accurately estimated from simulation. Second, the method typically results in fewer and better scaled parameters that can be more readily estimated from experiments. Finally, the simulation time can be significantly reduced without sacrificing the accuracy of the solution.

Stochastic kinetic modeling of vesicular stomatitis virus intracellular growth.

By building kinetic models of biological networks one may advance the development of new modeling approaches while gaining insights into the biology. We focus here on building a stochastic kinetic model for the intracellular growth of vesicular stomatitis virus (VSV), a well-studied virus that encodes five genes. The essential network of VSV reactions creates challenges to stochastic simulation owing to (i) delayed reactions associated with transcription and genome replication, (ii) production of large numbers of intermediate proteins by translation, and (iii) the presence of highly reactive intermediates that rapidly fluctuate in their intracellular levels. We address these issues by developing a hybrid implementation of the model that combines a delayed stochastic simulation algorithm (DSSA) with Langevin equations to simulate the reactions that produce species in high numbers. Further, we employ a quasi-steady-state approximation (QSSA) to overcome the computational burden of small time steps caused by highly reactive species. The simulation is able to capture experimentally observed patterns of viral gene expression. Moreover, the simulation suggests that early levels of a low-abundance species, VSV L mRNA, play a key role in determining the production level of VSV genomes, transcripts, and proteins within an infected cell. Ultimately, these results suggest that stochastic gene expression contribute to the distribution of virus progeny yields from infected cells.

Implications of decoupling the intracellular and extracellular levels in multi-level models of virus growth.

Virus infections are characterized by two distinct levels of detail: the intracellular level describing how viruses hijack the host machinery to replicate, and the extracellular level describing how populations of virus and host cells interact. Deterministic, population balance models for viral infections permit incorporation of both the intracellular and extracellular levels of information. In this work, we identify assumptions that lead to exact, selective decoupling of the interaction between the intracellular and extracellular levels, effectively permitting solution of first the intracellular level, and subsequently the extracellular level. This decoupling leads to (1) intracellular and extracellular models of viral infections that have been previously reported and (2) a significant reduction in the computational expense required to solve the model. However, the decoupling restricts the behaviors that can be modeled. Simulation of a previously reported multi-level model demonstrates this decomposition when the intracellular level of description consists of numerous reaction events. Additionally, examples demonstrate that viruses can persist even when the intracellular level of description cannot sustain a steady-state production of virus (i.e., has only a trivial equilibrium). We expect the combination of this modeling framework with experimental data to result in a quantitative, systems-level understanding of viral infections and cellular antiviral strategies that will facilitate controlling both these infections and antiviral strategies.

Image-guided modeling of virus growth and spread.

Although many tools of cellular and molecular biology have been used to characterize single intracellular cycles of virus growth, few culture methods exist to study the dynamics of spatially spreading viruses over multiple generations. We have previously developed a method that addresses this need by tracking the spread of focal infections using immunocytochemical labeling and digital imaging. Here, we build reaction-diffusion models to account for spatio-temporal patterns formed by the spreading viral infection front as well as data from a single cycle of virus growth (one-step growth). Systems with and without the interferon-mediated antiviral response of the host cells are considered. Dynamic images of the spreading infections guide iterative model refinement steps that lead to reproduction of all of the salient features contained in the images, not just the velocity of the infection front. The optimal fits provide estimates for key parameters such as virus-host binding and the production rate of interferon. For the examined data, highly-lumped infection models that ignore the one-step growth dynamics provide a comparable fit to models that more accurately account for these dynamics, highlighting the fact that increased model complexity does not necessarily translate to improved fit. This work demonstrates how model building can facilitate the interpretation of experiments by highlighting contributions from both biological and methodological factors.

Comparison of Parameter Estimation Methods in Stochastic Chemical Kinetic Models: Examples in Systems Biology.

Stochastic chemical kinetics has become a staple for mechanistically modeling various phenomena in systems biology. These models, even more so than their deterministic counterparts, pose a challenging problem in the estimation of kinetic parameters from experimental data. As a result of the inherent randomness involved in stochastic chemical kinetic models, the estimation methods tend to be statistical in nature. Three classes of estimation methods are implemented and compared in this paper. The first is the exact method, which uses the continuous-time Markov chain representation of stochastic chemical kinetics and is tractable only for a very restricted class of problems. The next class of methods is based on Markov chain Monte Carlo (MCMC) techniques. The third method, termed conditional density importance sampling (CDIS), is a new method introduced in this paper. The use of these methods is demonstrated on two examples taken from systems biology, one of which is a new model of single-cell viral infection. The applicability, strengths and weaknesses of the three classes of estimation methods are discussed. Using simulated data for the two examples, some guidelines are provided on experimental design to obtain more information from a limited number of measurements.

Two classes of quasi-steady-state model reductions for stochastic kinetics.

The quasi-steady-state approximation (QSSA) is a model reduction technique used to remove highly reactive species from deterministic models of reaction mechanisms. In many reaction networks the highly reactive intermediates (QSSA species) have populations small enough to require a stochastic representation. In this work we apply singular perturbation analysis to remove the QSSA species from the chemical master equation for two classes of problems. The first class occurs in reaction networks where all the species have small populations and the QSSA species sample zero the majority of the time. The perturbation analysis provides a reduced master equation in which the highly reactive species can sample only zero, and are effectively removed from the model. The reduced master equation can be sampled with the Gillespie algorithm. This first stochastic QSSA reduction is applied to several example reaction mechanisms (including Michaelis-Menten kinetics) [Biochem. Z. 49, 333 (1913)]. A general framework for applying the first QSSA reduction technique to new reaction mechanisms is derived. The second class of QSSA model reductions is derived for reaction networks where non-QSSA species have large populations and QSSA species numbers are small and stochastic. We derive this second QSSA reduction from a combination of singular perturbation analysis and the Omega expansion. In some cases the reduced mechanisms and reaction rates from these two stochastic QSSA models and the classical deterministic QSSA reduction are equivalent; however, this is not usually the case.

State estimation approach for determining composition and growth rate of Si1-xGex chemical vapor deposition utilizing real-time ellipsometric measurements.

We present an algorithm that simultaneously deduces from real-time ellipsometric measurements both the growth rate and the composition of Si1-xGex films deposited via chemical vapor deposition. The heart of the algorithm is a dynamic, first-principles model of the deposition system and the ellipsometric sensor. The model predicts the ellipsometric parameters psi and Delta during film growth. An extended Kalman filter is developed that utilizes the sensor model and infers both the growth rate and the Ge composition of the deposited film in real time. Two simulations demonstrating the effectiveness of the algorithm are evaluated.

Stochastic simulation of catalytic surface reactions in the fast diffusion limit.

The master equation of a lattice gas reaction tracks the probability of visiting all spatial configurations. The large number of unique spatial configurations on a lattice renders master equation simulations infeasible for even small lattices. In this work, a reduced master equation is derived for the probability distribution of the coverages in the infinite diffusion limit. This derivation justifies the widely used assumption that the adlayer is in equilibrium for the current coverages and temperature when all reactants are highly mobile. Given the reduced master equation, two novel and efficient simulation methods of lattice gas reactions in the infinite diffusion limit are derived. The first method involves solving the reduced master equation directly for small lattices, which is intractable in configuration space. The second method involves reducing the master equation further in the large lattice limit to a set of differential equations that tracks only the species coverages. Solution of the reduced master equation and differential equations requires information that can be obtained through short, diffusion-only kinetic Monte Carlo simulation runs at each coverage. These simulations need to be run only once because the data can be stored and used for simulations with any set of kinetic parameters, gas-phase concentrations, and initial conditions. An idealized CO oxidation reaction mechanism with strong lateral interactions is used as an example system for demonstrating the reduced master equation and deterministic simulation techniques.

On the origins of approximations for stochastic chemical kinetics.

This paper considers the derivation of approximations for stochastic chemical kinetics governed by the discrete master equation. Here, the concepts of (1) partitioning on the basis of fast and slow reactions as opposed to fast and slow species and (2) conditional probability densities are used to derive approximate, partitioned master equations, which are Markovian in nature, from the original master equation. Under different conditions dictated by relaxation time arguments, such approximations give rise to both the equilibrium and hybrid (deterministic or Langevin equations coupled with discrete stochastic simulation) approximations previously reported. In addition, the derivation points out several weaknesses in previous justifications of both the hybrid and equilibrium systems and demonstrates the connection between the original and approximate master equations. Two simple examples illustrate situations in which these two approximate methods are applicable and demonstrate the two methods' efficiencies.