Rate kernel theory for pseudo-first-order kinetics of diffusion-influenced reactions and application to fluorescence quenching kinetics

Theoretical foundation of rate kernel equation approaches for diffusion-influenced chemical reactions is presented and applied to explain the kinetics of fluorescence quenching reactions. A many-body master equation is constructed by introducing stochastic terms, which characterize the rates of chemical reactions, into the many-body Smoluchowski equation. A Langevin-type of memory equation for the density fields of reactants evolving under the influence of time-independent perturbation is derived. This equation should be useful in predicting the time evolution of reactant concentrations approaching the steady state attained by the perturbation as well as the steady-state concentrations. The dynamics of fluctuation occurring in equilibrium state can be predicted by the memory equation by turning the perturbation off and consequently may be useful in obtaining the linear response to a time-dependent perturbation. It is found that unimolecular decay processes including the time-independent perturbation can be incorporated into bimolecular reaction kinetics as a Laplace transform variable. As a result, a theory for bimolecular reactions along with the unimolecular process turned off is sufficient to predict overall reaction kinetics including the effects of unimolecular reactions and perturbation. As the present formulation is applied to steady-state kinetics of fluorescence quenching reactions, the exact relation between fluorophore concentrations and the intensity of excitation light is derived.     

Oscillations in non-mass action kinetics models of biochemical reaction networks arising from pairs of subnetworks

It is well known that oscillations in models of biochemical reaction networks can arise as a result of a single negative cycle. On the other hand, methods for finding general network conditions for potential oscillations in large biochemical reaction networks containing many cycles are not well developed. A biochemical reaction network with any number of species is represented by a simple digraph and is modeled by an ordinary differential equation (ODE) system with non-mass action kinetics. The obtained graph-theoretic condition generalizes the negative cycle condition for oscillations in ODE models to the existence of a pair of subnetworks, where each subnetwork contains an even number of positive cycles. The technique is illustrated with a model of genetic regulation.     

Approximation of slow attracting manifolds in chemical kinetics by trajectory-based optimization approaches

Many common kinetic model reduction approaches are explicitly based on inherent multiple time scales and often assume and directly exploit a clear time scale separation into fast and slow reaction processes. They approximate the system dynamics with a dimension-reduced model after eliminating the fast modes by enslaving them to the slow ones. The corresponding restrictive assumption of full relaxation of fast modes often renders the resulting approximation of slow attracting manifolds inaccurate as a representation of the reduced model and makes the numerical solution of the nonlinear "reduction equations" particularly difficult in many cases where the gap in intrinsic time scales is not large enough. We demonstrate that trajectory optimization approaches can avoid such severe restrictions by computing numerical solutions that correspond to "maximally relaxed" dynamical modes in a suitable sense. We present a framework of trajectory-based optimization for model reduction in chemical kinetics and a general class of reduction criteria characterizing the relaxation of chemical forces along reaction trajectories. These criteria can be motivated geometrically exploiting ideas from differential geometry and fundamental physics and turn out to be highly successful in example applications. Within this framework, we provide results for the computational approximation of slow attracting low-dimensional manifolds in terms of families of optimal trajectories for a six-component hydrogen combustion mechanism.     

Multiple‐Well,multiple‐path unimolecular reaction systems. I. MultiWell computer program suite

Unimolecular reaction systems in which multiple isomers undergo simultaneous reactions via multiple decomposition reactions and multiple isomerization reactions are of fundamental interest in chemical kinetics. The computer program suite described here can be used to treat such coupled systems, including the effects of collisional energy transfer (weak collisions). The program suite consists of MultiWell, which solves the internal energy master equation for complex unimolecular reactions systems; DenSum, which calculates sums and densities of states by an exact‐count method; MomInert, which calculates external principal moments of inertia and internal rotation reduced moments of inertia; and Thermo, which calculates equilibrium constants and other thermodynamics quantities. MultiWell utilizes a hybrid master equation approach, which performs like an energy‐grained master equation at low energies and a continuum master equation in the vibrational quasicontinuum. An adaptation of Gillespie's exact stochastic method is used for the solution. The codes are designed for ease of use. Details are presented of various methods for treating weak collisions with virtually any desired collision step‐size distribution and for utilizing RRKM theory for specific unimolecular rate constants. © 2001 John Wiley & Sons, Inc. Int J Chem Kinet 33: 232–245, 2001     

The generalized transition state method

Current theories of unimolecular reaction rates are based on the transition state method which replaces internal reactant dynamics by an assumption of internal equilibrium. The present work is devoted to the development of generalized transition state method which allows effects such as nonergodicity and non-exponential decay to be accounted for within a simple theoretical framework. The derivation is quantum mechanical and not limited by any weak perturbation assumption. An effective hamiltonian is constructed for the reactant dynamics. The loss of amplitude due to reaction is accounted for by a dissipative term in the hamiltonian which is obtained on a phenomenological basis. The diagonalization of the hamiltonian allows the decay of reactant state to be predicted. The decay information is then used to set up a non-markovian master equation which in turn yields the rate coefficient for the reaction. The accuracy of the method is tested in one-dimensional model calculations in which particular attention is paid to decay by quantum mechanical tunneling through a potential barrier.     

The functional equation truncation method for approximating slow invariant manifolds: a rapid method for computing intrinsic low-dimensional manifolds

A slow manifold is a low-dimensional invariant manifold to which trajectories nearby are rapidly attracted on the way to the equilibrium point. The exact computation of the slow manifold simplifies the model without sacrificing accuracy on the slow time scales of the system. The Maas-Pope intrinsic low-dimensional manifold (ILDM) [Combust. Flame 88, 239 (1992)] is frequently used as an approximation to the slow manifold. This approximation is based on a linearized analysis of the differential equations and thus neglects curvature. We present here an efficient way to calculate an approximation equivalent to the ILDM. Our method, called functional equation truncation (FET), first develops a hierarchy of functional equations involving higher derivatives which can then be truncated at second-derivative terms to explicitly neglect the curvature. We prove that the ILDM and FET-approximated (FETA) manifolds are identical for the one-dimensional slow manifold of any planar system. In higher-dimensional spaces, the ILDM and FETA manifolds agree to numerical accuracy almost everywhere. Solution of the FET equations is, however, expected to generally be faster than the ILDM method.     

Reduction of chemical reaction networks using quasi-integrals

We present a numerical method to identify possible candidates for quasi-stationary manifolds in complex reaction networks governed by systems of ordinary differential equations. Inspired by singular perturbation theory, we examine the ratios of certain components of the reaction rate vector. Those ratios that rapidly approach a nearly constant value define a slow manifold for the original flow in terms of quasi-integrals, that is, functions that are nearly constant along the trajectories. The dimensionality of the original system is thus effectively reduced without reliance on a priori knowledge of the different time scales in the system. We also demonstrate the relation of our approach to singular perturbation theory which, in its simplest form, is just the well-known quasi-steady-state approximation. In two case studies, we apply our method to oscillatory chemical systems: the 6-dimensional hemin-hydrogen peroxide-sulfite pH oscillator and a 10-dimensional mechanistic model for the peroxidase-oxidase (PO) reaction system. We conjecture that the presented method is especially suited for a straightforward reduction of higher dimensional dynamical systems where analytical methods fail to identify the different time scales associated with the slow invariant manifolds present in the system.     

Relaxation of chemical potential and a generalized diffusion equation

The effect of slow structural relaxation in a solvent of high viscosity on the chemical potential driving the diffusion of penetrant molecules is described by a generalized diffusion equation with a memory term. The linearized version of this equation is solved for some special cases, and the correlation function of concentration fluctuations in thermodynamic equilibrium is calculated. As a result of the memory term, for very slow relaxation two different stages of the diffusion process can be distinguished.     

The invariant constrained equilibrium edge preimage curve method for the dimension reduction of chemical kinetics

This work addresses the construction and use of low-dimensional invariant manifolds to simplify complex chemical kinetics. Typically, chemical kinetic systems have a wide range of time scales. As a consequence, reaction trajectories rapidly approach a hierarchy of attracting manifolds of decreasing dimension in the full composition space. In previous research, several different methods have been proposed to identify these low-dimensional attracting manifolds. Here we propose a new method based on an invariant constrained equilibrium edge (ICE) manifold. This manifold (of dimension nr) is generated by the reaction trajectories emanating from its (nr-1)-dimensional edge, on which the composition is in a constrained equilibrium state. A reasonable choice of the nr represented variables (e.g., nr "major" species) ensures that there exists a unique point on the ICE manifold corresponding to each realizable value of the represented variables. The process of identifying this point is referred to as species reconstruction. A second contribution of this work is a local method of species reconstruction, called ICE-PIC, which is based on the ICE manifold and uses preimage curves (PICs). The ICE-PIC method is local in the sense that species reconstruction can be performed without generating the whole of the manifold (or a significant portion thereof). The ICE-PIC method is the first approach that locally determines points on a low-dimensional invariant manifold, and its application to high-dimensional chemical systems is straightforward. The "inputs" to the method are the detailed kinetic mechanism and the chosen reduced representation (e.g., some major species). The ICE-PIC method is illustrated and demonstrated using an idealized H2O system with six chemical species. It is then tested and compared to three other dimension-reduction methods for the test case of a one-dimensional premixed laminar flame of stoichiometric hydrogen/air, which is described by a detailed mechanism containing nine species and 21 reactions. It is shown that the error incurred by the ICE-PIC method with four represented species is small across the whole flame, even in the low temperature region.     

Exact solutions for the entropy production rate of several irreversible processes

We investigate thermal conduction described by Newton's law of cooling and by Fourier's transport equation and chemical reactions based on mass action kinetics where we detail a simple example of a reaction mechanism with one intermediate. In these cases we derive exact expressions for the entropy production rate and its differential. We show that at a stationary state the entropy production rate is an extremum if and only if the stationary state is a state of thermodynamic equilibrium. These results are exact and independent of any expansions of the entropy production rate. In the case of thermal conduction we compare our exact approach with the conventional approach based on the expansion of the entropy production rate near equilibrium. If we expand the entropy production rate in a series and keep terms up to the third order in the deviation variables and then differentiate, we find out that the entropy production rate is not an extremum at a nonequilibrium steady state. If there is a strict proportionality between fluxes and forces, then the entropy production rate is an extremum at the stationary state even if the stationary state is far away from equilibrium.     

Computing minimal entropy production trajectories: an approach to model reduction in chemical kinetics

Advanced experimental techniques in chemistry and physics provide increasing access to detailed deterministic mass action models for chemical reaction kinetics. Especially in complex technical or biochemical systems the huge amount of species and reaction pathways involved in a detailed modeling approach call for efficient methods of model reduction. These should be automatic and based on a firm mathematical analysis of the ordinary differential equations underlying the chemical kinetics in deterministic models. A main purpose of model reduction is to enable accurate numerical simulations of even high dimensional and spatially extended reaction systems. The latter include physical transport mechanisms and are modeled by partial differential equations. Their numerical solution for hundreds or thousands of species within a reasonable time will exceed computer capacities available now and in a foreseeable future. The central idea of model reduction is to replace the high dimensional dynamics by a low dimensional approximation with an appropriate degree of accuracy. Here I present a global approach to model reduction based on the concept of minimal entropy production and its numerical implementation. For given values of a single species concentration in a chemical system all other species concentrations are computed under the assumption that the system is as close as possible to its attractor, the thermodynamic equilibrium, in the sense that all modes of thermodynamic forces are maximally relaxed except the one, which drives the remaining system dynamics. This relaxation is expressed in terms of minimal entropy production for single reaction steps along phase space trajectories.     

A probability generating function method for stochastic reaction networks

In this paper we present a probability generating function (PGF) approach for analyzing stochastic reaction networks. The master equation of the network can be converted to a partial differential equation for PGF. Using power series expansion of PGF and Padé approximation, we develop numerical schemes for finding probability distributions as well as first and second moments. We show numerical accuracy of the method by simulating chemical reaction examples such as a binding-unbinding reaction, an enzyme-substrate model, Goldbeter-Koshland ultrasensitive switch model, and G(2)/M transition model.     

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.     

Analysis of the approximate slow invariant manifold method for reactive flow equations

The Approximate Slow Invariant Manifold method of Singh, Powers and Paolucci is a useful method for addressing model reduction in systems of reactive flow equations. It exploits separations of time scales between slow and fast species, and it generalizes the Intrinsic Low-Dimensional Manifold method, which was developed for model reduction in the context of reaction kinetics, to systems with diffusive and active transport. In this article, we present a mathematical analysis of the Approximate Slow Invariant Manifold method in the context of systems of reaction–diffusion equations with slow and fast reaction kinetics. Beginning with systems of two species (one slow and one fast), and then treating general systems with multiple slow and fast species, we explicitly determine the accuracy of the Approximate Slow Invariant Manifold method. We find that it is correct up to and including the terms of first order in the small parameter that measures the separation of the kinetics time scales, and that it captures many of the terms at second order, as well. Our analysis includes precise statements of the errors at second order, and we find that these are proportional to the slow components of the reaction–diffusion equation, as well as to the curvature of the critical manifold. We illustrate the results analytically on two prototypical examples, the Michaelis–Menten–Henri model with diffusion of the slow species and the Davis–Skodje model in which both the slow and fast species diffuse.     

Chemical kinetics with reagent dispersion in single-line flow-injection systems

Simultaneous chemical kinetics and dispersion kinetics in single-line flow-injection manifold can be modeled by using total rate laws for dispersion and chemical reaction. The model combines a (non-equilibrium) mass-transfer rate coefficient with residence-time theory based on extended tanks-in-series. Reagent dispersion is found to play a key role in defining a zone of mixing. Double product peaks and narrow analyte peaks observed in the literature are predicted by the model. When applied to a slow chemical reaction under conditions of constant length of mixing, predicted and experimentally observed response curves agree well.     

Modeling of nonlinear boundary value problems in enzyme-catalyzed reaction diffusion processes

A mathematical model of steady state mono-layer potentiometric biosensor is developed. The model is based on non stationary diffusion equations containing a non linear term related to Michaelis-Menten kinetics of the enzymatic reaction. This paper presents a complex numerical method (He’s variational iteration method) to solve the non-linear differential equations that describe the diffusion coupled with a Michaelis-Menten kinetics law. Approximate analytical expressions for substrate concentration and corresponding current response have been derived for all values of saturation parameter α and reaction diffusion parameter K using variational iteration method. These results are compared with available limiting case results and are found to be in good agreement. The obtained results are valid for the whole solution domain.     

Reducing a chemical master equation by invariant manifold methods

We study methods for reducing chemical master equations using the Michaelis-Menten mechanism as an example. The master equation consists of a set of linear ordinary differential equations whose variables are probabilities that the realizable states exist. For a master equation with s(0) initial substrate molecules and e(0) initial enzyme molecules, the manifold can be parametrized by s(0) of the probability variables. Fraser's functional iteration method is found to be difficult to use for master equations of high dimension. Building on the insights gained from Fraser's method, techniques are developed to produce s(0)-dimensional manifolds of larger systems directly from the eigenvectors. We also develop a simple, but surprisingly effective way to generate initial conditions for the reduced models.