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

Published on

04-Sep-2016View

212Download

0

Embed Size (px)

Transcript

<ul><li><p>Stiff stochastic differential equations which a positive denite Jacobian is maintained to deal with possible stiffness, is proposed</p><p>s) arling pn theh CLEropen</p><p>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.</p><p>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</p><p>0021-9991/$ - see front matter 2011 Elsevier Inc. All rights reserved.</p><p> Corresponding author.E-mail addresses: saswatid@rishi.serc.iisc.ernet.in (S. Dana), raha@serc.iisc.ernet.in (S. Raha).</p><p>Journal of Computational Physics 230 (2011) 88138834</p><p>Contents lists available at SciVerse ScienceDirect</p><p>Journal of Computational Physicsdoi:10.1016/j.jcp.2011.07.032sive 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.</p><p>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).</p><p>1.1. Computational issues</p><p>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</p><p>1. Introduction</p><p>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.</p><p> 2011 Elsevier Inc. All rights reserved.</p><p>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</p><p>Saswati Dana, Soumyendu Raha Supercomputer Education and Research Centre, Indian Institute of Science, Bangalore 560012, India</p><p>a r t i c l e i n f o</p><p>Article history:Received 3 December 2010Received in revised form 17 June 2011Accepted 26 July 2011Available online 10 September 2011</p><p>Keywords:Chemical Langevin equations</p><p>a b s t r a c t</p><p>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</p><p>journal homepage: www.elsevier .com/locate / jcp</p></li><li><p>cost ationalmatedintegra reaseventsmalltation</p><p>Ho</p><p>t</p><p>adaptcoefc</p><p>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><p>aiXt</p><p>p. In more computationally efcient CLE mod-</p><p>els [19], the m noise channels are reduced to a Wiener process with dimension equal to the number of independent reaction</p><p>directions, say, l 6 m and B is obtained as an ll matrix constructed, for example, in a mTJVdiagaXtVT</p><p>qform as in Con-</p><p>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:</p><p>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:</p><p>LetVTdiagaUZtV</p><p>q: GZt : Rl ! Rll:</p><p>Let W be l-dimensional Wiener process: Then</p><p>dZt RaZtdt RGZtdWt;</p><p>3</p><p>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</p><p>BBT mTdiagaXtm 2j1</p><p>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</p><p>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</p><p>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.</p><p>2. The chemical Langevin equation models</p><p>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</p><p>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</p></li><li><p>2.1. No</p><p>Let (Then</p><p>The prdiffusically</p><p>j(j)</p><p>j k Xt ; if each Xt P 0;</p><p>holds f~</p><p>so tha</p><p>S. Dana, S. Raha / Journal of Computational Physics 230 (2011) 88138834 8815EVt20;tk 1 kX0k22 p=2</p><p>Z t0ELVs ds 6 1 kX0k22</p><p> p=2 2p max</p><p>j2 1;...; ~mf gjjkmk2F</p><p>Z t0EVs ds 9i1 j1 i1 j1 j2f1;...; ~mg</p><p>t for a stopping time tk, 0 < tk 6 T, the variation of the Lyapunov function, LV, satises6 p 1 kXsk22 p=21 Xl</p><p>XisX~m</p><p>mj;iajXs Xl X~m</p><p>mj;i 2</p><p>ajXs 6 p 1 kXsk22 p=2</p><p>2 max jjkmk2F ; 8 pp 22</p><p>1 kXsk22 p=22X</p><p>i1aiXs</p><p>Xj1</p><p>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.</p><p>Proof. Set the Lyapunov function Vt 1 kXtk22p=2 for t 2 [0,T], T 0; j 2 f1; . . . ;mg; 7</p><p>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</p><p>In the present work we estimate the stiffness of a CLE system as the local exponential growth or decay of the states over anite time interval. This is physically consistent with dening the stiffness of CLE (1) as a sharp change in the state of a spe-cies over a small enough but non-zero time interval.</p><p>2.2.1. Exponential growthIn the following we estimate the local exponential growth (e.g., see [18]) bound for a CLE system.</p><p>Lemma 1. Let pP 2 and X0 2 V be A0-measurable and is known almost surely (i.e., with probability 1), whereV fx 2 Rl : 1 > xi P 0 8i 2 f1; . . . ; lgg. Suppose the linear growth bounda Xt i2Rj0; otherwise:</p><p>>: 6are reactants in the jth reaction channel havingmj reactants, k is the rate constant for the jth reaction channel and wi: > 0 isthe number of moles reacting of the i.th species, then</p><p>j Q i wi i8> xi > 0 8i 2 f1; . . . ; lgg. Consistent with the shutting down of a reac-tion channel when any one of its reactants has become depleted (reached zero state value), and with the chemical kinetics,a(Xt) is assumed to be in the following form. If Rj fi1; . . . ; im g where i. 2 {1, . . . , l} is the set of indices of the species which0 if Xit 6 0:</p><p>the CLEs (1) can be written as</p><p>Xt I X0 Z t0mTaXsds</p><p>Z t0BXsdWs</p><p> : 5</p><p>opensities are usually in the form of the product of the states of the species and from (1) and (3) we can see that theion coefcients in CLEs involve square root of linear combinations of the propensities. In view of this it will not be phys-inconsistent to assume that the components of the propensities vector a(Xt) and those of the diffusion coefcient B(Xt)</p><p>lIiXt Xit if X</p><p>it > 0; 4n-negativity and propensities</p><p>IXt : Rl ! Rl dened as</p></li><li><p>by Itis esti</p><p>p=2 p=2</p><p>2.2.2.</p><p>ma 1, we discuss the effect of this phenomenon on stiffness of a CLE system.Let</p><p>cientlchannels active) and in [s,T], T 0 leading to large maxj jj. Similar to (10) or (11), if we derive the growth of EkXtk22, the local exponen-non-negativity constraint (1b).Moreover, due to the square root diffusion coefcient in the CLEs, the local linear growth bound of the propensities inhannel has the largest linear growth rate constant then, by Lemma 1 the upper bounded exponential growth rate forE system changes from pmaxj2f1;...;m1gjjkmk2F to pjmkmk2F . This illustrates the relationship between stiffness and theby Its formula and assumption (7), as in Lemma 1. Applying Gronwall inequality it can be readily seen that if themth reac-that the mth reaction channel switches on at s having been inactive until then. Then, we compute</p><p>EVt20;s 1 kX0k22 p=2</p><p>Z t0ELVs ds 6 1 kX0k22</p><p> p=2 2p max</p><p>j2f1;...;m1gjjkmk2F</p><p>Z s0EVs ds</p><p>and</p><p>EVt2s;T E 1 kXsk22 p=2</p><p>Z t</p><p>ELVs ds 6 E 1 kXsk22 p=2</p><p> 2p max jjkmkj2FZ T</p><p>EVs dss be a nite non-zero stopping time such that in [0,s), m 1 reaction channels are active in a system (i.e., l is suf-y large so that when one component of Xt is zero for t 2 [0,s), the other non-depleted species keep the m 1 reactionNon-negativity of the CLEs physically imply that if a species is depleted its state becomes zero and a reaction channel inwhich it is a reactant shuts down, until the same species is produced in some other active reaction channel. In view of Lem-Switching of reaction channels and stiffnessE 1 kXtk22 6 1 kX0k22 exp 2p maxj2f1;...; ~mg</p><p>jjkmk2F t : 10</p><p>for t 2 (0,T]. Hence the result. h</p><p>Physically, if the reaction system consists of reaction channels with widely scattered rate constants (speeds), then afterscaling the CLE (1) with a suitable reaction rate constant, we will have large (scaled) rate constants corresponding to thefaster reactions. The scaled propensities of the relatively faster reaction channels have a larger linear growth bound, i.e., alarger jwhich depends on the rate constant. Applying (10) to the scaled CLE, we see that a larger maxj2f1;...; ~mgjj leads to shar-per exponential growth, i.e., stiffness.</p><p>Remark 2. Also, one can observe that under the contractivity condition</p><p>xTmTa 12tracemTdiagam < 0</p><p>in place of (7), an exponential decay bound in the form of</p><p>E 1 kXtk22 </p><p>6 1 kX0k22 </p><p>exp ~jt ; 11</p><p>where ~j is a suitable constant dependent on the reaction rate constants, can be obtained, proceeding as in Lemma 1.s formula. Then, by Gronwall inequality and letting the stopping time sequence tk? T, k?1, the exponential growthmated as </p></li><li><p>tn1 s j j </p><p>1 m l ll</p><p>j1 ;j2 j1 ;j2</p><p>n n1 n n n1 n</p><p>where</p><p>m m ! ! !</p><p>ations, Pn is calculated such that the real parts of the eigenvalues of</p><p>are st</p><p>where</p><p>S. Dana, S. Raha / Journal of Computational Physics 230 (2011) 88138834 8817provides a low pass lter mechanism giving the method its large mean-square stability region and in turn, its effectivenesson stiff SDEs.I bcaPn @</p><p>@xfh</p><p>Xmj11</p><p>Xmj21</p><p>Dcj1 ;j2Ij1 ;j2n</p><p> ! !</p><p>rictly positive. The algorithmic coefcients are computed as</p><p>c 2q 1</p><p>12; 18a</p><p>b 1q 12 ; 18b</p><p>Y0 h2 @f@x</p><p>f </p><p>X0X0: h2K0 0; 17b</p><p>Jui;j : uieXn1 uiXn1eXjn1 Xjn1 for any u : Rl ! Rl; j~Xn1 Xnj > 0;kPnk2 1 almost surely; 17aXn1 Xn f n bc aPnD/n 12 bc</p><p> Yn; 16a</p><p>Yn1 1cD/n 11c</p><p> Yn; 16b</p><p>1 a > 1 is a user-selectable real parameter; Pn is a real matrix P evaluated at the nth time step such that,FIS-a method [22] is given asn n1 n</p><p>D/n : DfnhXmj11</p><p>Xmj21</p><p>Dcj1 ;j2n Ij1 ;j2n 15fn : f Xn; Bn : BXn; cn : c Xn;Df : f f ; DB : B B ; Dcj1 ;j2 : cj1 ;j2 cj1 ;j2;B : fb ; . . . ; b g : R ! R ;ri : @</p><p>@xiith component of the gradient operator;</p><p>Lj : bjTr; cj1 ;j2 : Lj1bj2 : Rl ! Rl; 13f n : f Xn : f Xnh BXnDWn </p><p>Xmj11</p><p>Xmj21</p><p>cj1 ;j2nXnIj1 ;j2;...</p></li></ul>