Physically consistent simulation of mesoscale chemical kinetics: The non-negative FIS-α method

  • Published on

  • View

  • Download

Embed Size (px)


  • Stiff stochastic differential equations which a positive denite Jacobian is maintained to deal with possible stiffness, is proposed

    s) arling pn theh CLEropen

    produce non-negative sample paths than a lower strong order convergent method. Besides, a higher order method is usefulwhen the non-global Lipschitz continuity of the diffusion coefcients can lead to loss of the order of accuracy of thenumerical approximation.

    CLE models themselves being approximation of the CMEs matching only in the rst two moments [8,19], numericalschemes with a sufcient order of weak convergence of the rst two moments to the discretized CLE is usually adequate

    0021-9991/$ - see front matter 2011 Elsevier Inc. All rights reserved.

    Corresponding author.E-mail addresses: (S. Dana), (S. Raha).

    Journal of Computational Physics 230 (2011) 88138834

    Contents lists available at SciVerse ScienceDirect

    Journal of Computational Physicsdoi:10.1016/ as the integrator may be forced to take very small time steps. Additionally, the processes in the CLE systems are non-negative ones. The presence of square root of the propensities in the diffusion coefcients of the CLEs is a source of numericaldifculty with the Jacobians of implicit methods, particularly at time points where a reactant species approaches depletion.

    The aim of the present work is to numerically compute individual non-negative sample paths consistent with the CLEsthat are It stochastic differential equations (SDEs).

    1.1. Computational issues

    When a CLE system has non-negative solutions over the simulation interval over which the reaction channels switch offand on, a higher strong order convergent numerical scheme is more likely (considering Doobs martingale inequality) toMeso-scale kineticsImplicit method of numerical integration

    1. Introduction

    Chemical Langevin equations (CLEing some of the intra-cellular signalsizes of various species participate imaster equation (CME) [21], of whicreaction rate constant values in the pand analysed. The method is illustrated with the computation of active Protein Kinase Cresponse in the Protein Kinase C pathway.

    2011 Elsevier Inc. All rights reserved.

    e used to model mesoscale reaction kinetics in many bio-chemical systems includ-athways [17] where a medium sized population [8] of molecules of comparablechemical reactions. Stiffness (sharp changes in short time scales) in the chemicals are approximations, is seen in the solution of CLE systems having wide ranges ofsities. Numerical integration of such CLEs systems can be computationally expen-Physically consistent simulation of mesoscale chemical kinetics: Thenon-negative FIS-a method

    Saswati Dana, Soumyendu Raha Supercomputer Education and Research Centre, Indian Institute of Science, Bangalore 560012, India

    a r t i c l e i n f o

    Article history:Received 3 December 2010Received in revised form 17 June 2011Accepted 26 July 2011Available online 10 September 2011

    Keywords:Chemical Langevin equations

    a b s t r a c t

    Biochemical pathways involving chemical kinetics in medium concentrations (i.e., at meso-scale) of the reacting molecules can be approximated as chemical Langevin equations (CLE)systems. We address the physically consistent non-negative simulation of the CLE samplepaths as well as the issue of non-Lipschitz diffusion coefcients when a species approachesdepletion and any stiffness due to faster reactions. The non-negative Fully Implicit Stochas-tic a (FIS a) method in which stopped reaction channels due to depleted reactants aredeleted until a reactant concentration rises again, for non-negativity preservation and in

    journal homepage: www.elsevier .com/locate / jcp

  • cost ationalmatedintegra reaseventsmalltation




    8814 S. Dana, S. Raha / Journal of Computational Physics 230 (2011) 88138834and B mTdiag aXtp in the standard CLE as derived in [8]; where the entry (i, i) of the diagonal matrix diag aXtp isgiven by the square root of the ith component of propensity vector, i.e.,


    p. In more computationally efcient CLE mod-

    els [19], the m noise channels are reduced to a Wiener process with dimension equal to the number of independent reaction

    directions, say, l 6 m and B is obtained as an ll matrix constructed, for example, in a mTJVdiagaXtVT

    qform as in Con-

    struction 4 of [19]. For computation, the CLE can then be further reduced in state space by a suitable linear transformation(cf. Theorem 3.6 in [19]). All these alternative formulations of the standard CLEs constructed by noise channel and state spacereductions are equivalent and the species states obtained from them match in all the moments with the standard CLE. Ingeneral it may be seen that the following construction gives an irreducible CLE model:

    Let lP mP l rankm:Let mT URVT be the thin singular value decomposition of mT:Zt 2 Rl; Xt UZt ; Zt UTXt ; VTaUZt : aZt : Rl ! Rl:


    q: GZt : Rl ! Rll:

    Let W be l-dimensional Wiener process: Then

    dZt RaZtdt RGZtdWt;


    can be easily veried to be the most reduced (in both state space and noise channels) equivalent of the standard CLE and is inthe form of (1a), although the computational use of a reduced CLE as per the constructions in [19] or as (3) may depend onthe computational cost of the linear transformations and/or of a sufciently accurate thin singular value decomposition and/or of the square root of the positive denite matrices. Overall (1) is sufciently general to include both standard and reducedforms of the CLEs and our numerical method is analysed based on the same.ed to ltration fAt ; t P 0g on a common probability space X;A;P. The initial data X0 is A0-measurable. The diffusionient B satises

    BBT mTdiagaXtm 2j1

    Xit P 0; 8i 2 f1; . . . ; lg; 1bwhere X is a non-negative real valued stochastic process, Xt 2 Rl is the number (approximated as real) of molecules or theconcentration of a species at time t 2 T : 0; T R; a : Rl ! Rl is the vector of the propensity functions and a is positiveincreasing over Rl; b

    j is the jth column of the diffusion coefcient B : Rl ! Rlm; m 2 Rml (usually having integer entriesand, in particular, is {1,0,1}ml for bi-molecular reactions) is the stoichiometric matrix relating ith reaction channel(row) to jth species state (column), and W = {W , tP 0} is an m-dimensional Wiener process with independent channels,dXt mTaXtdt Xm

    bjXtdW jt ; 1atency of a sample path, such as the non-negativity and oscillation periods of species concentrations in each realization, isof interest to the present study. Hence, we explore a numerical method nominally convergent with strong order 1.0 for SDEswith globally Lipschitz coefcients, as a trade-off between the strong order 1/2 methods and the computational cost of meth-ods with strong order greater than 1.0.

    2. The chemical Langevin equation models

    We consider the CLE system represented as an It stochastic differential equation (SDE) system with a non-negativityconstraint on the state of the species participating in the chemical reaction network:nd the error in mean hitting time of the CLE approximations with respect to the CMEs. The error is inversely propor-to the reaction rate, i.e., mean hitting times of faster reactions which can make a CME simulation costly is approxi-better (cf. Theorem 6.1 in [24]). A numerical method with higher strong order convergence may be desirable for

    ating the CLE approximation of the said class of bio-chemical networks in order to meet computational objectives overonably moderate number of realizations and to reduce the error in estimation of time intervals between biochemicals (e.g., oscillations and depletion of concentrations). For example, in [17,23] Milstein schemes have been used with verytime steps and the results have been compared with exact CME simulations for trade-off between error and compu-al effort.wever higher order methods involve the increased computational cost of multiple It integrals. The physical consis-for computation. However in some applications the mean time to hitting a particular level of a concentration may be of inter-est (e.g., in [17] the active Protein Kinase C oscillation is studied). For large biochemical pathway reaction networks involvinga wide range of reaction rates, the exact stochastic simulation of CMEs are computationally very expensive [21]. In somesuch systems, even the implicit tau-leaping schemes (e.g., the method described in [21]) may prove expensive. These sys-tems, such as the pathway in [17], are simulated using the CLE approximations as a trade-off between the computational

  • 2.1. No

    Let (Then

    The prdiffusically


    j k Xt ; if each Xt P 0;

    holds f~

    so tha

    S. Dana, S. Raha / Journal of Computational Physics 230 (2011) 88138834 8815EVt20;tk 1 kX0k22 p=2

    Z t0ELVs ds 6 1 kX0k22

    p=2 2p max

    j2 1;...; ~mf gjjkmk2F

    Z t0EVs ds 9i1 j1 i1 j1 j2f1;...; ~mg

    t for a stopping time tk, 0 < tk 6 T, the variation of the Lyapunov function, LV, satises6 p 1 kXsk22 p=21 Xl


    mj;iajXs Xl X~m

    mj;i 2

    ajXs 6 p 1 kXsk22 p=2

    2 max jjkmk2F ; 8 pp 22

    1 kXsk22 p=22X



    mi;jXjs ! !species in the CLE over a nite non-zero time interval is proportional tomaxj2f1;...; ~mgjj, where m is number of active reaction chan-nels over the given time interval.

    Proof. Set the Lyapunov function Vt 1 kXtk22p=2 for t 2 [0,T], T 0; j 2 f1; . . . ;mg; 7

    or all x 2 V. Let Xt 2 V for all t 2 [0,T] for the CLE (1). Then the upper bound of the exponential growth rate of the states of the2.2. Stiffness of the CLEs