Muruhan Rathinam

Stiffness in stochastic chemically reacting systems: The implicit tau-leaping method

We show how stiffness manifests itself in the simulation of chemical reactions at both the continuous-deterministic level and the discrete-stochastic level. Existing discrete stochastic simulation methods, such as the stochastic simulation algorithm and the (explicit) tau-leaping method, are both exceedingly slow for such systems. We propose an implicit tau-leaping method that can take much larger time steps for many of these problems.

A New Look at Proper Orthogonal Decomposition

We investigate some basic properties of the proper orthogonal decomposition (POD) method as it is applied to data compression and model reduction of finite dimensional nonlinear systems. First we provide an analysis of the errors involved in solving a nonlinear ODE initial value problem using a POD reduced order model. Then we study the effects of small perturbations in the ensemble of data from which the POD reduced order model is constructed on the reduced order model. We explain why in some applications this sensitivity is a concern while in others it is not. We also provide an analysis of computational complexity of solving an ODE initial value problem and study the computational savings obtained by using a POD reduced order model. We provide several examples to illustrate our theoretical results.

Configuration Flatness of Lagrangian Systems Underactuated by One Control

Lagrangian control systems that are differentially flat with flat outputs that depend only on configuration variables are said to be configuration flat. We provide a complete characterization of configuration flatness for systems with n degrees of freedom and n-1 controls whose range of control forces only depends on configuration and whose Lagrangian has the form of kinetic energy minus potential. The method presented allows us to determine if such a system is configuration flat and, if so, provides a constructive method for finding all possible configuration flat outputs. Our characterization relates configuration flatness to Riemannian geometry. We illustrate the method with two examples.

Differential flatness and absolute equivalence

In this paper we give a formulation of differential flatness-a concept originally introduced by Fliess, Levine, Martin, and Rouchon (1992)-in terms of absolute equivalence between exterior differential systems. Systems which are differentially flat have several useful properties which can be exploited to generate effective control strategies for nonlinear systems. The original definition of flatness was given in the context of differential algebra, and required that all mappings be meromorphic functions. Our formulation of flatness does not require any algebraic structure and allows one to use tools from exterior differential systems to help characterize differentially flat systems. In particular, we show that in the case of single input control systems (i.e., codimension 2 Pfaffian systems), a system is differentially flat if and only if it is feedback linearizable via static state feedback. However, in higher codimensions feedback linearizability and flatness are not equivalent: one must be careful with the role of time as well the use of prolongations which may not be realizable as dynamic feedbacks in a control setting. Applications of differential flatness to nonlinear control systems and open questions are also discussed.<<ETX>>

Efficient computation of parameter sensitivities of discrete stochastic chemical reaction networks

Parametric sensitivity of biochemical networks is an indispensable tool for studying system robustness properties, estimating network parameters, and identifying targets for drug therapy. For discrete stochastic representations of biochemical networks where Monte Carlo methods are commonly used, sensitivity analysis can be particularly challenging, as accurate finite difference computations of sensitivity require a large number of simulations for both nominal and perturbed values of the parameters. In this paper we introduce the common random number (CRN) method in conjunction with Gillespies stochastic simulation algorithm, which exploits positive correlations obtained by using CRNs for nominal and perturbed parameters. We also propose a new method called the common reaction path (CRP) method, which uses CRNs together with the random time change representation of discrete state Markov processes due to Kurtz to estimate the sensitivity via a finite difference approximation applied to coupled reaction paths that emerge naturally in this representation. While both methods reduce the variance of the estimator significantly compared to independent random number finite difference implementations, numerical evidence suggests that the CRP method achieves a greater variance reduction. We also provide some theoretical basis for the superior performance of CRP. The improved accuracy of these methods allows for much more efficient sensitivity estimation. In two example systems reported in this work, speedup factors greater than 300 and 10,000 are demonstrated.

Differential Flatness and Absolute Equivalence of Nonlinear Control Systems

This paper presents a formulation of differential flatness---a concept originally introduced by Fliess, Levine, Martin, and Rouchon---in terms of absolute equivalence between exterior differential systems. Systems that are differentially flat have several useful properties that can be exploited to generate effective control strategies for nonlinear systems. The original definition of flatness was given in the context of differential algebra and required that all mappings be meromorphic functions. The formulation of flatness presented here does not require any algebraic structure and allows one to use tools from exterior differential systems to help characterize differentially flat systems. In particular, it is shown that, under regularity assumptions and in the case of single input control systems (i.e., codimension 2 Pfaffian systems), a system is differentially flat if and only if it is feedback linearizable via static state feedback. In higher codimensions our approach does not allow one to prove that feedback linearizability about an equilibrium point and flatness are equivalent: one must be careful with the role of time as well as the use of prolongations that may not be realizable as dynamic feedback in a control setting. Applications of differential flatness to nonlinear control systems and open questions are also discussed.

Consistency and Stability of Tau-Leaping Schemes for Chemical Reaction Systems

We develop a theory of local errors for the explicit and implicit tau-leaping methods for simulating stochastic chemical systems, and we prove that these methods are first-order consistent. Our theory provides local error formulae that could serve as the basis for future stepsize control techniques. We prove that, for the special case of systems with linear propensity functions, both tau-leaping methods are first-order convergent in all moments. We provide a stiff stability analysis of the mean of both leaping methods, and we confirm that the implicit method is unconditionally stable in the mean for stable systems. Finally, we give some theoretical and numerical examples to illustrate these results.

The numerical stability of leaping methods for stochastic simulation of chemically reacting systems

Tau-leaping methods have recently been proposed for the acceleration of discrete stochastic simulation of chemically reacting systems. This paper considers the numerical stability of these methods. The concept of stochastic absolute stability is defined, discussed, and applied to the following leaping methods: the explicit tau, implicit tau, and trapezoidal tau.

Dynamic Iteration Using Reduced Order Models: A Method for Simulation of Large Scale Modular Systems

We describe a new iterative method, dynamic iteration using reduced order models (DIRM), for simulation of large scale modular systems using reduced order models that preserve the interconnection structure. This method may be compared to the waveform relaxation technique; however, unlike DIRM, waveform relaxation does not take advantage of model reduction techniques. The DIRM method involves simulating in turn each subsystem connected to model reduced versions of the other subsystems. The data from this simulation is then used to update the reduced model for that particular subsystem. We provide analytical results on convergence and accuracy of the DIRM method as well as numerical examples that demonstrate the success of DIRM and verify the analysis.

A pathwise derivative approach to the computation of parameter sensitivities in discrete stochastic chemical systems

Characterizing the sensitivity to infinitesimally small perturbations in parameters is a powerful tool for the analysis, modeling, and design of chemical reaction networks. Sensitivity analysis of networks modeled using stochastic chemical kinetics, in which a probabilistic description is used to characterize the inherent randomness of the system, is commonly performed using Monte Carlo methods. Monte Carlo methods require large numbers of stochastic simulations in order to generate accurate statistics, which is usually computationally demanding or in some cases altogether impractical due to the overwhelming computational cost. In this work, we address this problem by presenting the regularized pathwise derivative method for efficient sensitivity analysis. By considering a regularized sensitivity problem and using the random time change description for Markov processes, we are able to construct a sensitivity estimator based on pathwise differentiation (also known as infinitesimal perturbation analysis) that is valid for many problems in stochastic chemical kinetics. The theoretical justification for the method is discussed, and a numerical algorithm is provided to permit straightforward implementation of the method. We show using numerical examples that the new regularized pathwise derivative method (1) is able to accurately estimate the sensitivities for many realistic problems and path functionals, and (2) in many cases outperforms alternative sensitivity methods, including the Girsanov likelihood ratio estimator and common reaction path finite difference method. In fact, we observe that the variance reduction using the regularized pathwise derivative method can be as large as ten orders of magnitude in certain cases, permitting much more efficient sensitivity analysis than is possible using other methods.