A partial-propensity formulation of the stochastic simulation algorithm for chemical reaction networks with delays

Several real-world systems, such as gene expression networks in biological cells, contain coupled chemical reactions with a time delay between reaction initiation and completion. The non-Markovian kinetics of such reaction networks can be exactly simulated using the delay stochastic simulation algorithm (dSSA). The computational cost of dSSA scales with the total number of reactions in the network. We reduce this cost to scale at most with the smaller number of species by using the concept of partial reaction propensities. The resulting delay partial-propensity direct method (dPDM) is an exact dSSA formulation for well-stirred systems of coupled chemical reactions with delays. We detail dPDM and present a theoretical analysis of its computational cost. Furthermore, we demonstrate the implications of the theoretical cost analysis in two prototypical benchmark applications. The dPDM formulation is shown to be particularly efficient for strongly coupled reaction networks, where the number of reactions is much larger than the number of species.     

Extracting the mechanisms and kinetic models of complex reactions from atomistic simulation data

Determining reaction mechanisms and kinetic models, which can be used for chemical reaction engineering and design, from atomistic simulation is highly challenging. In this study, we develop a novel methodology to solve this problem. Our approach has three components: (1) a procedure for precisely identifying chemical species and elementary reactions and statistically calculating the reaction rate constants; (2) a reduction method to simplify the complex reaction network into a skeletal network which can be used directly for kinetic modeling; and (3) a deterministic method for validating the derived full and skeletal kinetic models. The methodology is demonstrated by analyzing simulation data of hydrogen combustion. The full reaction network comprises 69 species and 256 reactions, which is reduced into a skeletal network of 9 species and 30 reactions. The kinetic models of both the full and skeletal networks represent the simulation data well. In addition, the essential elementary reactions and their rate constants agree favorably with those obtained experimentally. © 2019 Wiley Periodicals, Inc.     

Multinomial tau-leaping method for stochastic kinetic simulations

We introduce the multinomial tau-leaping (MtauL) method for general reaction networks with multichannel reactant dependencies. The MtauL method is an extension of the binomial tau-leaping method where efficiency is improved in several ways. First, tau-leaping steps are determined simply and efficiently using a priori information and Poisson distribution-based estimates of expectation values for reaction numbers over a tentative tau-leaping step. Second, networks are partitioned into closed groups of reactions and corresponding reactants in which no group reactant set is found in any other group. Third, product formation is factored into upper-bound estimation of the number of times a particular reaction occurs. Together, these features allow larger time steps where the numbers of reactions occurring simultaneously in a multichannel manner are estimated accurately using a multinomial distribution. Furthermore, we develop a simple procedure that places a specific upper bound on the total reaction number to ensure non-negativity of species populations over a single multiple-reaction step. Using two disparate test case problems involving cellular processes--epidermal growth factor receptor signaling and a lactose operon model--we show that the tau-leaping based methods such as the MtauL algorithm can significantly reduce the number of simulation steps thus increasing the numerical efficiency over the exact stochastic simulation algorithm by orders of magnitude.     

Molecular device design based on chemical reaction networks: state feedback controller,static pre-filter,addition gate control system and full-dimensional state observer

Journal of Mathematical Chemistry - For the modeling and operation of biological computing, chemical reaction networks (CRNs) constitute an ideal programming paradigm for the simulation of various...     

Preprogramming Complex Hydrogel Responses using Enzymatic Reaction Networks

The creation of adaptive matter is heavily inspired by biological systems. However, it remains challenging to design complex material responses that are governed by reaction networks, which lie at the heart of cellular complexity. The main reason for this slow progress is the lack of a general strategy to integrate reaction networks with materials. Herein we use a systematic approach to preprogram the response of a hydrogel to a trigger, in this case the enzyme trypsin, which activates a reaction network embedded within the hydrogel. A full characterization of all the kinetic rate constants in the system enabled the construction of a computational model, which predicted different hydrogel responses depending on the input concentration of the trigger. The results of the simulation are in good agreement with experimental findings. Our methodology can be used to design new, adaptive materials of which the properties are governed by reaction networks of arbitrary complexity.     

Efficient computation of transient solutions of the chemical master equation based on uniformization and quasi-Monte Carlo

A quasi-Monte Carlo method for the simulation of discrete time Markov chains is applied to the simulation of biochemical reaction networks. The continuous process is formulated as a discrete chain subordinate to a Poisson process using the method of uniformization. It is shown that a substantial reduction of the number of trajectories that is required for an accurate estimation of the probability density functions (PDFs) can be achieved with this technique. The method is applied to the simulation of two model problems. Although the technique employed here does not address the typical stiffness of biochemical reaction networks, it is useful when computing the PDF by replication. The method can also be used in conjuncture with hybrid methods that reduce the stiffness.     

An electrochemical adaptation of Ellman's test

The electrochemically initiated reaction of thiols with N,N-diethyl-p-phenylenediamine has been coupled with an existing colorimetric sensing reaction developed by Ellman as a means of providing an electrochemical adaptation of the latter whereby the total thiol species present in a sample can be determined. The detection methodology has been proven to be robust with a linear range for cysteine from 2-120 microM, a limit of detection of 1.17 microM and has shown selectivity against a wide range of potential interferences. The efficiency of the methodology has been examined in the determination and recovery of thiol species in three growth tissue media, which contain a number of common biological interferences.     

Determination of complex reaction mechanisms. Analysis of chemical, biological and genetic networks

We present several methods of determining, not guessing, complex chemical reaction mechanisms and their functions. One method is based on the theory of correlation functions of measured time series of concentrations of chemical species; another is on measurements of temporal responses of concentrations to various perturbations of arbitrary magnitude; a third deals with the analysis of oscillatory systems; a fourth is on the use of genetic algorithms to determine functions of chemical reaction networks. All methods are applicable to chemical, biochemical, and biological reaction systems and to genetic networks and systems biology. The methods depend on the design of appropriate experiments on the whole system and corresponding theories for interpretation that lead to information on the causal chemical connectivity of species, on reaction pathways, on reaction mechanisms, on control centers in the system, and on functions of the system. The first three methods require no assumption of a model or hypothesis, nor extensive calculations, unlike the interpretation of measurements made on a gene network at only one time.     

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.     

Modeling signal transduction networks: a comparison of two stochastic kinetic simulation algorithms

Computational efficiency of stochastic kinetic algorithms depend on factors such as the overall species population, the total number of reactions, and the average number of nodal interactions or connectivity in a network. These size measures of the network model can have a significant impact on computational efficiency. In this study, two scalable biological networks are used to compare the size scaling efficiencies of two popular and conceptually distinct stochastic kinetic simulation algorithms--the random substrate method of Firth and Bray (FB), and the Gillespie algorithm as implemented using the Gibson-Bruck method (GGB). The arithmetic computational efficiencies of these two algorithms, respectively, scale with the square of the total species population and the logarithm of the total number of active reactions. The two scalable models considered are the size scalable model (SSM), a four compartment reaction model for a signal transduction network involving receptors with single phosphorylation binding sites, and the variable connectivity model (VCM), a single compartment model where receptors possess multiple phosphorylation binding sites. The SSM has fixed species connectivity while the connectivity between species in VCM increases with the number of phosphorylation sites. For SSM, we find that, as the total species population is increased over four orders of magnitude, the GGB algorithm performs significantly better than FB for all three SSM compartment models considered. In contrast, for VCM, we find that as the overall species population decreases while the number of phosphorylation sites increases (implying an increase in network linkage) there exists a crossover point where the computational demands of the GGB method exceed that of the FB.     

A reduction method for multiple time scale stochastic reaction networks

In this paper we develop a reduction method for multiple time scale stochastic reaction networks. When the transition-rate matrix between different states of the species is available, we obtain systems of reduced equations, whose solutions can successively approximate, to any degree of accuracy, the exact probability that the reaction system be in any particular state. For the case when the transition-rate matrix is not available, one needs to rely on the chemical master equation. For this case, we obtain a corresponding reduced master equation with first-order accuracy. We illustrate the accuracy and efficiency of both approaches by simulating several motivating examples and comparing the results of our simulations with the results obtained by the exact method. Our examples include both linear and nonlinear reaction networks as well as a three time scale stochastic reaction-diffusion model arising from gene expression.     

Redox active metal-induced oxidative stress in biological systems

A number of studies performed on biological systems have shown that redox-active metals such as iron and copper as well as other transition metals can undergo redox cycling reactions and produce reactive free radicals termed also reactive oxygen species (ROS) or reactive nitrogen species (RNS). The most representative examples of ROS and RNS are the superoxide anion radical and nitric oxide, respectively, both playing a dual role in biological systems. At low/moderate concentrations of ROS and RNS, they can be involved in many physiological roles such as defense against infectious agents, involvement in a number of cellular signaling pathways and other important biological processes. On the other hand, at high concentrations, ROS and RNS can be important mediators of damage to biomolecules involving DNA, membrane lipids, and proteins. One of the most damaging ROS occurring in biological systems is the hydroxyl radical formed via the decomposition of hydrogen peroxide catalyzed by traces of iron, copper and other metals (the Fenton reaction). The hydroxyl radical is known to react with the DNA molecule, forming 8-OH-Guanine adduct, which is a good biomarker of oxidative stress of an organism and a potential biomarker of carcinogenesis. This review discusses the role of iron and copper in uncontrolled formation of ROS leading to various human diseases such as cancer, cardiovascular disease, and neurological disorders (Alzheimer’s disease and Parkinson’s disease). A discussion is devoted to the various protective antioxidant networks against the deleterious action of free radicals. Metal-chelation therapy, which is a modern pharmacotherapy used to chelate redox-active metals and remove toxic metals from living systems to avoid metal poisoning, is also discussed.     

Rate constant estimation for C1 to C4 alkyl and alkoxyl radical decomposition

Rate coefficients for alkyl and alkoxy radical decomposition are important in combustion, biological, and atmospheric processes. In this paper, rate constant expressions for C₁? C₄ alkyl and alkoxy radicals decomposition via β‐scission are recommended based on the reverse, exothermic reaction, the addition of a hydrogen atom or an alkyl radical to an olefin or carbonyl species with the decomposition reaction calculated using microscopic reversibility. The rate expressions have been estimated based on a wide‐range study of available experimental data. Rate coefficients for hydrogen atom and alkyl radical addition to an olefin show a strong temperature curvature. In addition, it is found that there is a correlation between the activation energy for addition and (i) the type of atom undergoing addition and (ii) whether this radical adds to the internal or terminal carbon atom of the olefin. Rate coefficients for alkoxy radical decomposition show a strong correlation to the ionization potential of the alkyl radical leaving group and on the enthalpy of reaction. It is shown that the activation energy for alkyl radical addition to a carbonyl species can be estimated as a function of the alkyl radical ionization potential and enthalpy of reaction. © 2006 Wiley Periodicals, Inc. Int J Chem Kinet 38: 250–275, 2006     

ReaxFF MD局部区域反应追踪与物理性质可视化分析

Recently, the application of ReaxFF based reactive molecular dynamics simulation (ReaxFF MD) in complex processes of pyrolysis, oxidation and catalysis has attracted considerable attention. The analysis of the simulation results of these processes is challenging owing to the complex chemical reactions involved, coupled with their dynamic physical properties. VARxMD is a leading tool for the chemical reaction analysis and visualization of ReaxFF MD simulations, which allows the automated analysis of reaction sites to get overall reaction lists, evolution trends of reactants and products, and reaction networks of specified reactants and products. The visualization of the reaction details and dynamic evolution profiles are readily available for each reactant and product. Additionally, the detailed reaction sites of bond breaking and formation are available in 2D chemical structure diagrams and 3D structure views; for specified reactions, they are categorized on the basis of the chemical structures of the bonding sites or function groups in the reacting species. However, the current VARxMD code mainly focuses on global chemical reaction information in the simulation system of the ReaxFF MD, and is incapable of locally tracking the chemical reaction and physical properties in a 3D picked zone. This work extends the VARxMD from global analysis to a focused 3D zone picked interactively from the 3D visualization modules of VARxMD, as well as physical property analysis to complement reaction analysis. The analysis of reactions and physical properties can be implemented in three steps: picking and drawing a 3D zone, identifying molecules in the picked zone, and analyzing the reactions and physical properties of the picked molecules. A 3D zone can be picked by specifying the geometric parameters or drawing on a screen using a mouse. The picking of a cuboid or sphere was implemented using the VTK 3D view libraries by specifying geometric parameters. The interactive 3D zone picking was implemented using a combination of observer and command patterns in the VTK visualization paradigm. The chemical reaction tracking and dynamic radial distribution function (RDF) of the 3D picked zone was efficiently implemented by inheriting data obtained from the global analysis of VARxMD. The reaction tracking between coal particles in coal pyrolysis simulation and dynamic structure characterization of carbon rich cluster formation in the thermal decomposition of an energetic material are presented as application examples. The obtained detailed reactions between the coal particles and comparison of the reaction between the locally and globally picked areas in the cuboid are helpful in understanding the role of micro pores in coal particles. The carbon to carbon RDF analysis and comparison of the spherical region picked for the layered molecular clusters in the pyrolysis system of the TNT crystal model with the standard RDF of the 5-layer graphene demonstrate the extended VARxMD as a chemical structure characteristic tool for detecting the dynamic formation profile of carbon rich clusters in the pyrolysis of TNT. The extended capability of VARxMD for a 3D picked zone of a ReaxFF MD simulation system can be useful for interfacial reaction analysis in a catalysis system, hot spot formation analysis in the detonation of energetic material systems, and particularly the pyrolysis or oxidation processes of coal, biomass, polymers, hydrocarbon fuels, and energetic materials.     

Investigations of Pd-catalyzed ArX coupling reactions informed by reaction progress kinetic analysis

This Perspective highlights how the methodology of reaction progress kinetic analysis can provide a rapid and comprehensive kinetic profile of complex catalytic reaction networks under synthetically relevant conditions in a fraction of the number of experiments required by classical kinetic analysis. This approach relies on graphical manipulation of the extensive data sets available from accurate in situ monitoring of reaction progress under conditions where two concentration variables are changing simultaneously. A series of examples from Pd-catalyzed coupling reactions of aryl halides demonstrates how a wealth of kinetic information may be extracted from just three experiments in each case. Even before proposing a reaction mechanism, we can determine reaction orders in substrates, propose a resting state for the catalyst, and probe catalyst stability. Carrying out this kinetic analysis at the outset of a mechanistic investigation provides a framework for further work aimed at seeking a molecular-level understanding of the nature of the species within the catalytic cycle. To be considered plausible, any independent mechanistic proposal must be shown to be consistent with this global kinetic analysis.     

An accelerated algorithm for discrete stochastic simulation of reaction-diffusion systems using gradient-based diffusion and tau-leaping

Stochastic simulation of reaction-diffusion systems enables the investigation of stochastic events arising from the small numbers and heterogeneous distribution of molecular species in biological cells. Stochastic variations in intracellular microdomains and in diffusional gradients play a significant part in the spatiotemporal activity and behavior of cells. Although an exact stochastic simulation that simulates every individual reaction and diffusion event gives a most accurate trajectory of the system's state over time, it can be too slow for many practical applications. We present an accelerated algorithm for discrete stochastic simulation of reaction-diffusion systems designed to improve the speed of simulation by reducing the number of time-steps required to complete a simulation run. This method is unique in that it employs two strategies that have not been incorporated in existing spatial stochastic simulation algorithms. First, diffusive transfers between neighboring subvolumes are based on concentration gradients. This treatment necessitates sampling of only the net or observed diffusion events from higher to lower concentration gradients rather than sampling all diffusion events regardless of local concentration gradients. Second, we extend the non-negative Poisson tau-leaping method that was originally developed for speeding up nonspatial or homogeneous stochastic simulation algorithms. This method calculates each leap time in a unified step for both reaction and diffusion processes while satisfying the leap condition that the propensities do not change appreciably during the leap and ensuring that leaping does not cause molecular populations to become negative. Numerical results are presented that illustrate the improvement in simulation speed achieved by incorporating these two new strategies.     

A numerical characteristic method for probability generating functions on stochastic first-order reaction networks

We propose an efficient and accurate numerical scheme for solving probability generating functions arising in stochastic models of general first-order reaction networks by using the characteristic curves. A partial differential equation derived by a probability generating function is the transport equation with variable coefficients. We apply the idea of characteristics for the estimation of statistical measures, consisting of the mean, variance, and marginal probability. Estimation accuracy is obtained by the Newton formulas for the finite difference and time accuracy is obtained by applying the fourth order Runge–Kutta scheme for the characteristic curve and the Simpson method for the integration on the curve. We apply our proposed method to motivating biological examples and show the accuracy by comparing simulation results from the characteristic method with those from the stochastic simulation algorithm.     

A "partitioned leaping" approach for multiscale modeling of chemical reaction dynamics

We present a novel multiscale simulation approach for modeling stochasticity in chemical reaction networks. The approach seamlessly integrates exact-stochastic and "leaping" methodologies into a single partitioned leaping algorithmic framework. The technique correctly accounts for stochastic noise at significantly reduced computational cost, requires the definition of only three model-independent parameters, and is particularly well suited for simulating systems containing widely disparate species populations. We present the theoretical foundations of partitioned leaping, discuss various options for its practical implementation, and demonstrate the utility of the method via illustrative examples.     

An equation-free probabilistic steady-state approximation: dynamic application to the stochastic simulation of biochemical reaction networks

Stochastic chemical kinetics more accurately describes the dynamics of "small" chemical systems, such as biological cells. Many real systems contain dynamical stiffness, which causes the exact stochastic simulation algorithm or other kinetic Monte Carlo methods to spend the majority of their time executing frequently occurring reaction events. Previous methods have successfully applied a type of probabilistic steady-state approximation by deriving an evolution equation, such as the chemical master equation, for the relaxed fast dynamics and using the solution of that equation to determine the slow dynamics. However, because the solution of the chemical master equation is limited to small, carefully selected, or linear reaction networks, an alternate equation-free method would be highly useful. We present a probabilistic steady-state approximation that separates the time scales of an arbitrary reaction network, detects the convergence of a marginal distribution to a quasi-steady-state, directly samples the underlying distribution, and uses those samples to accurately predict the state of the system, including the effects of the slow dynamics, at future times. The numerical method produces an accurate solution of both the fast and slow reaction dynamics while, for stiff systems, reducing the computational time by orders of magnitude. The developed theory makes no approximations on the shape or form of the underlying steady-state distribution and only assumes that it is ergodic. We demonstrate the accuracy and efficiency of the method using multiple interesting examples, including a highly nonlinear protein-protein interaction network. The developed theory may be applied to any type of kinetic Monte Carlo simulation to more efficiently simulate dynamically stiff systems, including existing exact, approximate, or hybrid stochastic simulation techniques.     

关于副本交换分子动力学模拟复杂化学反应的研究

本文研究用温度副本交换分子动力学(T-REMD)和哈密顿副本交换分子动力学(H-REMD)方法模拟复杂化学反应的问题。使用具有不同活化能和反应能的简单置换反应模型,我们检验了上述两种方法用来预测反应平衡产物的效率和应用范围。T-REMD方法对具有适度活化能(约< 20 kcal·mol^-1)或者反应能量(< 3 kcal·mol^-1)的放热反应是有效的。由于在相空间的不完整采样,对于同时具有高活化能和反应能量的反应其模拟效率有严重障碍,并且对于吸热反应问题更为显着。另一方面,H-REMD对一系列具有不同活化能的反应能的模型表现出色,与T-REMD相比,H-REMD可以使用更少的副本获得优异的结果。