Documents

Bioinformatics-2009-Surovtsova-2816-23.pdf

Categories
Published
of 8
All materials on our website are shared by users. If you have any questions about copyright issues, please report us to resolve them. We are always happy to assist you.
Related Documents
Share
Description
[11:05 22/10/2009 Bioinformatics-btp451.tex] Page: 2816 2816–2823 BIOINFORMATICS ORIGINAL PAPER Vol. 25 no. 21 2009, pages 2816–2823 doi:10.1093/bioinformatics/btp451 Systems biology Accessible methods for the dynamic time-scale decomposition of biochemical systems Irina Surovtsova 1,∗ , Natalia Simus 1 , Thomas Lorenz 2 , Artjom König 1 , Sven Sahle 1 and Ursula Kummer 1 1 Department of Modeling of Biological Processes and 2 Interdisciplinary Center for Scientific Computing, Univers
Transcript
   BIOINFORMATICS ORIGINAL PAPER  Vol. 25 no. 21 2009, pages 2816–2823doi:10.1093/bioinformatics/btp451 Systems biology   Accessible methods for the dynamic time-scale decomposition ofbiochemical systems Irina Surovtsova 1 , ∗ , Natalia Simus 1 , Thomas Lorenz 2 , Artjom König 1 , Sven Sahle 1 andUrsula Kummer 1 1 Department of Modeling of Biological Processes and  2 Interdisciplinary Center for Scientific Computing,University of Heidelberg, Im Neuenheimer Feld 267, 69120 Heidelberg, Germany Received on December 11, 2008; revised on May 22, 2009; accepted on June 29, 2009 Advance Access publication July 24, 2009 Associate Editor: Trey Ideker  ABSTRACTMotivation:  The growing complexity of biochemical models asksfor means to rationally dissect the networks into meaningful andrather independent subnetworks. Such foregoing should ensurean understanding of the system without any heuristics employed.Important for the success of such an approach is its accessibilityand the clarity of the presentation of the results. Results:  In order to achieve this goal, we developed a method whichis a modification of the classical approach of time-scale separation.This modified method as well as the more classical approachhave been implemented for time-dependent application within thewidely used software COPASI. The implementation includes differentpossibilities for the representation of the results including 3D-visualization.  Availability:  The methods are included in COPASI which is free foracademic use and available at www.copasi.org. Contact:  irina.surovtsova@bioquant.uni-heidelberg.de Supplementary information:  Supplementary data are available at Bioinformatics  online. 1 INTRODUCTION Animportantcomputationalaspectinsystemsbiologyisthefactthatthe increasing size and complexity of studied biochemical systemslead not only to experimental results which are hard to understand,but also to computational results which are not easy to comprehend.Here, the so-called complexity reduction aims for two differentdirections: increasing the speed of simulations and dissecting thebiochemical network into smaller subsystems that can be studiedindependently. The first one is usually achieved by reducing thenumber of equations mathematically. This has often been tackled byanalyzingthepresenceofawiderangeofcharacteristictimescalesinthe respective systems.Among the most prominent approaches builton this concept are the Computational Singular Perturbation (CSP)method of Lam and Goussis (1994) and the so-called Intrinsic low-dimensional manifolds (ILDM) method of Maas and Pope (1992).Several variants of these two methods have been used successfully,e.g. in atmospheric and combustion chemistry modeling (see e.g.Schmidt  et al . 1998). The CSP method was recently applied to the ∗ To whom correspondence should be addressed. model reduction of circardian rhythm in Drosophila by Goussisand Najm (2006). Besides, Lovrics  et al.  (2006) have used time-scale analysis to determine the stability of the budding yeast-cycletrajectories.A smaller number of equations, though, does not guaranteethe reduction of biochemical species in the system, since manydifferent species might contribute to one and the same transformedequation. Therefore, a likewise or even more important aspectis to dissect the biochemical system into several modules thatcan be studied independently. This is needed to understand theinterplay of specific subsystems. Most of the past approachesto complexity reduction of biochemical systems have focusedmainly on methods studying the steady-state behavior of the system(Kauffman  et al ., 2002; Price  et al ., 2003), or on dissecting thesystem based on its network topology using heuristic rules (Holme et al ., 2003; Schuster  et al ., 2002). The first approach is helpfulfor biochemical systems that indeed can be expected to displaysteady-state behavior. Nevertheless, most organisms will not displaysteady states, but rather transient behavior of different kinds atall times. Moreover, for large and complex biochemical reactionnetworks an additional requirement is that reduction methodscan be applied in an automatic way according to well-definedprinciples. For this purpose an automatic method based on a time-dependent ILDM was developed by Zobeley  et al.  (2005) forbiochemical systems.An additional analysis of species contributionto the slow dynamics has been introduced making the actualdecomposition of a system feasible (see Zobeley  et al. , 2005, aswell as Shaik   et al. , 2005).Another automatic complexity reductionmethod for biochemistry was presented recently by Lebiedz  et al. (2008). This combines dynamic sensitivity analysis, singular valuedecomposition and further computation of metabolite contributionto the active dynamics of the system.However, in addition to automation, it is equally important toensure accessibility of the respective methods, e.g. by integratingthemincommonlyavailableandfrequentlyusedsoftwarepackages.That goal is not trivial since the representation of the results, e.g.via visualization is quite difficult.This article presents the development of a new method, as wellas the implementation of this and an established method of Zobeley et al.  (2005) for the systematic and automated reduction proceduresthat may be applied to any explicit biochemical mechanism. 2816  © The Author 2009. Published by Oxford University Press. All rights reserved. For Permissions, please email: journals.permissions@oxfordjournals.org   b  y g u e  s  t   on O c  t   o b  e r 1  6  ,2  0 1 4 h  t   t   p :  /   /   b i   oi  nf   or m a  t  i   c  s  . oxf   or  d  j   o ur n a l   s  . or  g /  D o wnl   o a  d  e  d f  r  om    Dynamic time-scale decomposition of biochemical systems Both methods rely on the principles of the classical ILDMmethod which decomposes biochemical systems computationallywith respect to time scales. However, in comparison to the classicalILDM method applied, e.g. to chemical systems, the methods havean additional focus on the reduction of the underlying biochemicalnetwork in a time-dependent manner. The first approach has alreadybeen published in Zobeley  et al . (2005) and also used in Shaik   et al .(2005) and Surovtsova  et al . (2006). It is now implemented in auser-friendly way in the software framework COPASI (Hoops  et al .,2006). The implementation includes a visualization of the results inthe form of tables and 3D bar graphs.The second algorithm—also implemented in COPASI—usesessentially the same numerical tools as the classical ILDM.Its main goal, however, is to distinguish between ‘slow’ and‘fast’ metabolites—instead of ‘modes’. Hence, the final resultof the analysis can be formulated without linear transformationof the reaction system. Both methods can be applied toarbitrary biochemical reaction systems and work independently of assumptions about the specific dynamic regime.The article is organized in the following way: first, we describethe algorithms for time-scale separation as implemented in COPASI(Hoops  et al ., 2006). Then, we proceed to describe details of therespective implementation and finally, both methods are applied toseveral biochemical examples like calcium oscillations and yeastglycolysis. We also discuss the tolerance criterion used for detectingthe dimension of the fast space. 2 TIME-SCALE DECOMPOSITION: ALGORITHM Both algorithms presented here rely on the presence of a wide rangeof characteristic time scales in biological systems and are based onthe local analysis of the Jacobian, which is partitioned into fast andslow components at the initial point of a user-chosen interval. Weassume that the dynamics is determined by a system of   n  ordinarydifferentialequations(ODE)togetherwithagivenstate c 0  attime t  0 :dd t  c ( t  ) = f  ( c ( t  )) ,  c ( t  0 ) = c 0 ∈ R n (1)Moreover, the metabolite concentration  c , reaction rates  v  andstoichiometry  N  are related by: d c / d t  = N · v . We assume that all linear dependencies due to conservationrelationships have already been removed from the system. 2.1 Distinction between ‘fast’ and ‘slow’: QSSA In this article, the terms ‘fast’ and ‘slow’ are used for bothmetabolites [described by their time-dependent concentrations( c 1 ( t  ) ... c n ( t  )) ∈ R n ] and ‘modes’ [resulting from appropriatelyconstructed transformations of the metabolite concentrations as inDeuflhard and Heroth (1996) and Zobeley  et al . (2005)]. Theirbasic notion, however, is always closely related to the question howrapidly they change in time. The rate of change is usually describedby the time derivative.The starting point is partitioning the reaction system into slowand fast contributions. Such a splitting is related to a singularperturbation description of the transformed ODE system   dd t  x slow  =  g slow  x slow , x fast  ε ·  dd t  x fast  =  g fast  x slow , x fast  Here,  ε  is a singular perturbation parameter. The concentrations of the fast species change with time, but these species can be describedby algebraic relations instead of differential equations. The quasisteady-state assumption (QSSA) yields the associated differential–algebraic system (DAE) for the reduced problem by replacing ε ·  dd t   x fast ≈ 0 . The step to DAEs, however, provides difficulties with regard to theinitial conditions. In general, the given initial state ( x 0 , slow , x 0 , fast )does not satisfy the algebraic equation  g fast ( x 0 , slow , x 0 , fast ) = 0 . For this reason, we replace the known initial value  x 0 , fast  by anapproximation  x alg 0 , fast  consistent with the algebraic equation 0 = g fast ( x 0 , slow , x alg 0 , fast ) .  (2)These consequences of QSSAprovide our mathematical characteri-zation of ‘slow’and ‘fast’components: we regard those componentsas ‘fast’for which the QSSAcan be justified (on the basis of a fixederror tolerance).The other components are considered as ‘slow’and,they are still described by the full ODE. 2.2 Time-dependent classical ILDM Step 1:As a starting point of analysis, a linearization of Equation (1)with respect to the state vector  c 0  is performed:   dd t   c ( t  )  =  f   c 0  + J c 0 ·  c ( t  ) − c 0  c ( t  0 )  =  c 0  J c 0  Jacobian of   f   in  c 0  Step 2: An orthogonal similarity transformation is applied to theJacobian  J c 0 . The resulting matrix  S  has real Schur form, i.e. it is ablocks upper triangular matrix (Golub and van Loan, 1996): Q − 1 · J c 0 · Q = S =  S slow  S coup 0  S fast  The matrix  S  is reordered according to characteristic time-scales(TS):  τ  i =− (1 / Re λ i ) by means of Givens rotations ( λ i  denote theeigenvalues of the Jacobian).Step 3: Using the solution  Z  of the Sylvester equation S slow  Z − Z S fast =− S coup we realize that the transformed Jacobian is decoupled additionally: T − 1 · J c 0 · T =  S slow  00  S fast  where T = Q  Id n +   0  Z 0 0  ,  T − 1 =  Id n −   0  Z 0 0  Q − 1 Step 4: Applying  T − 1 to the state  c  and reaction rate  f   results in adecoupled representation of the system dynamics: x =  x slow x fast  = T − 1 · c ,  g =  g slow g fast  = T − 1 · f  ( T · )Step 5: The number  r   of ‘slow’ modes plays a crucial role in thisapproach. The criterion suggested by Deuflhard and Heroth (1996)is adapted (compare also Zobeley  et al . 2005): for a given errortolerancetol > 0,thenumber r   ofslowmodesisiterativelydecreased 2817   b  y g u e  s  t   on O c  t   o b  e r 1  6  ,2  0 1 4 h  t   t   p :  /   /   b i   oi  nf   or m a  t  i   c  s  . oxf   or  d  j   o ur n a l   s  . or  g /  D o wnl   o a  d  e  d f  r  om    I.Surovtsova et al. as long as the corresponding decomposition satisfies ε ·  g slow  x 0 , slow , x 0 , fast  − g slow  x 0 , slow , x alg 0 , fast  < tol (3)where the vector  x 0 , slow , x 0 , fast   is the transformed initial vector c 0  of ODE system (1) and  x 0 , slow , x alg 0 , fast   denotes its closestapproximation that is consistent with the algebraic Equation (2). ε = τ  r  + 1 = 1 / | Re λ r  + 1 | is a singular perturbation parameter.Step 6: The matrix  T  relates modes  x  to the srcinal variables  c .In order to complete our investigation, we perform the study of modes in terms of contributions of all species concentrations byanalyzing the entries of transformation matrices  T = ( β ik  ) 1 ≤ i , k  ≤ n and  T − 1 = ( γ  ik  ) 1 ≤ i , k  ≤ n  (see Surovtsova  et al . 2006; Zobeley  et al .2005).Vector ‘Slow space’ (‘Fast space’, respectively) indicates thecontribution of each concentration variable to the set of all slow(fast) modes:‘Slow space’( i ) :=  r k  = 1  γ  ik   n j , k  = 1  γ   jk  ,  ‘Fast space’( i ) :=  nk  = r  + 1  γ  ik   n j , k  = 1  γ   jk  2.3 Modified time-dependent ILDM method For proposing a selection of fast and slow metabolites, we useessentially the same numerical tools of local Schur decompositionof the Jacobian as for the classical variant above. The vector ‘Fastspace’ specifies the metabolites making the largest contributionsto the set of fast modes. Hence, it provides the next candidate of concentration which is then tested whether its additional nominationas fast metabolite still satisfies the criterion of error tolerance. Thealgorithm of our modified method could be described as follows:Step 1: linearize the dynamical system at  c 0 .Step 2: detect all time scales of the system and transform theODEs (sorted by absolute eigenvalues alias time scales) using Schurdecomposition.Step 3: analyze the Schur transformation matrix  Q  for obtaining thecontribution of each metabolite to the fast space. The vector ‘Fastspace’ specifies the metabolites making the largest contributions tothe set of fast modes.Step 4: check whether the QSSA for the most dominant metabolitestill satisfies the tolerance criterion, i.e. the dimension of fast space(and the number of fast metabolites) is increased as long as ε ·  f  slow  c 0 , slow , c 0 , fast  − f  slow  c 0 , slow , c alg 0 , fast  < tol (4)holds. Here,  c 0 , slow , c 0 , fast   denotes the initial state vector  c 0  of ODE (1) with respect to the current distinction between ‘slow’ and‘fast’ metabolites and, the vector  c 0 , slow , c alg 0 , fast   is its consistentinitial value of DAEs. As before,  ε = τ  r  + 1 = 1 / | Re λ r  + 1 | representsthe slowest time scale of the fast space. 2.3.1 Remark   In this modified ILDM method, we dispense withsolving the Sylvester equation and hence with decoupling  S slow  and S fast .Thusingeneral,wecannotexcludethattheslowmodesdependvery much on the exact values of fast modes. Nevertheless, the non-zero matrix  S coup  does not influence the number of fast metabolites.We also note that the non-zero matrix  S coup  does not influence thederivation of error estimate (see Deuflhard and Heroth, 1996). 2.4 Tolerance criterion for the new method Deuflhard and Heroth (1996) provide a mathematical justificationfor their tolerance criterion (3) on the basis of asymptotic analysis(with respect to the perturbation parameter  ε ).Since  ε  is not accessible in our case, we suggest an alternative justification of the tolerance criterion (4) whose mathematicalsimplicity might be regarded as an advantage. c 0 , slow + ( t  − t  0 ) · f  slow  c 0 , slow , c 0 , fast   is the (probably simplest)approximation of the exact concentration  c slow ( t  ) after a shortperiod  t  − t  0 ≥ 0 as the integral of the derivative easily reveals.Similarly,  c 0 , slow + ( t  − t  0 ) · f  slow  c 0 , slow , c alg 0 , fast   approximates thisconcentration resulting from the corresponding DAE and the initialstate  c 0 , slow , c alg 0 , fast  . Hence, the right-hand side of inequality (4)estimates the absolute error of the ‘slow’ metabolites up to time  ε. 2.5 Implementation in COPASI The Jacobian is calculated by using finite differences. The Schurtransformation and the solution of Sylvester equation are performedusing CLAPACK. The options for both methods are the same:ã Intervals: the user specifies the number of intervals to whichthe method is applied. COPASI performs the analysis at eachinitial point of these intervals. The interval size should be largein comparison with the time scale interesting for the user.ã Deuflhard tolerance: this parameter is a positive numericvalue specifying the maximal tolerated error of slow modes(Deuflhard and Heroth, 1996). The default value is 1 × 10 − 6 .We note that reasonable values of Deuflhard tolerance shouldbe of the same order as the time scale the user is interested in. 2.6 Results and graphical visualization in COPASI As described above results are displayed in two matrices and twovectors. The matrix ‘  Metabolites ’provides the contribution of eachmetabolite to every mode, whereas the table ‘  Modes ’ summarizesthe mode distribution for each metabolite.In the classical variant [called ILDM(Deuflhard)], the vector‘ Slow space ’ (‘ Fast space ’, respectively) indicates the relativecontribution of each concentration variable to the set of all slow(fast) modes. The metabolites with the largest contribution to thefast space could be supposed to be ‘fast’ and thus, its ODE can bereplaced by the corresponding algebraic equation (i.e. with the sameright-hand side). In the specific case that a subset of species does notcontribute to the slow space (but contributes only to the fast space),the time-scale decomposition results can be used for a dissection of the reaction network.The results of the modified variant consist of two similarmatrices and two similar vectors, which describe the contributionof metabolites to slow (fast) modes (spaces). However, thecorresponding numbers of slow and fast components now refer tothe distinction between slow and fast metabolites (not modes).All these matrices and vectors can either be exported to a text fileor displayed in the graphical user interface as tables. In this case, acolor coding is used where the numbers are additionally visualizedby different shades of color. This makes it easy to immediately spot,e.g. the most important contributions to a specific mode for a largemodel (where the result tables are correspondingly large).We also use 3D bar graphs for visualizing the matrices (see alsoFig. 5). These bar graphs can be turned and zoomed interactively. 2818   b  y g u e  s  t   on O c  t   o b  e r 1  6  ,2  0 1 4 h  t   t   p :  /   /   b i   oi  nf   or m a  t  i   c  s  . oxf   or  d  j   o ur n a l   s  . or  g /  D o wnl   o a  d  e  d f  r  om    Dynamic time-scale decomposition of biochemical systems Furthermore, single rows or columns of the matrix can behighlighted. An additional diagram shows the distribution of thetime scales of the different modes at chosen points of time (seeFig. 4).In the graphical user interface, it is very simple to switch betweenthe results for different points of time. Therefore, the user can easilyget an overview over the time-dependent changes of the time-scaleseparation. 2.7 ILDM and classical QSS and QE assumptions Two of the commonly used reduction methods, based on time-scale analysis, employ the quasi steady-state (QSS) approximationand quasi equilibrium (QE) assumption. QSSA identifies specieswhose production and destruction rate are in approximate balance,mathematically, it means that the right-hand side  f  i  of thecorresponding equation is zero. QE assumption corresponds toreactions whose forward and reverse rates are nearly equal.Nevertheless, the identification of such fast reactions and steady-state species is normally based on experience and intuition.The ILDM methods implemented in COPASI enable theautomatic decoupling of fast and slow processes. The analysis of transformation matrices allows the detection of QE reactions andQSS metabolites.The ( n − r  ) species dominating in the vector ‘Fast space’ inthe modified ILDM satisfy tolerance criterion (4). It justifiesthe QSS approximation for corresponding metabolites. The fastreactions could be discovered from the analysis of product betweenstoichiometric and transformation matrices in ‘classical ILDM’.The reactions which dominate the fast modes and which are lessimportant for the slow modes could be regarded as QE. 3 EXAMPLES 3.1 Michaelis–Menten kinetics We start our discussion with the simplest enzymatic reactionmechanism: Michaelis–Menten irreversible kinetics, described by Substrate Enzyme Complex Product Enzyme S   +  E   k  1 ⇐⇒ k  − 1 C   k  2 −→  P  +  E  PalssonandLightfoot(1984)exploredthepropertiesofthemodelvia improved scaling and linearization about a characteristic steadystate. Here, we use their analytical approach for testing our time-scale separation methods. A slightly different scaling of modelvariables is used: u = S S  0 ,  v = C  E  0 ,  e =  E  E  0 , τ  = k  1  E  0 t  , where  E  0  and  S  0  are the corresponding initial concentrations.Substitution of these scaled variables into the differential equationsleads to the following scaled equations:d u d τ  = − u ( τ  )  +  u ( τ  ) + Q 1 + S  t   v ( τ  )  dv d  τ  =  u ( τ  )  −  u ( τ  ) + Q  v ( τ  ) Table 1.  Michael–Menten model:   → 0Classical ILDM Modified methodTS C S C S4.2626 39.68 60.32 0.23 99.770.0034 99.77 0.23 99.77 0.23 The table displays metabolite contribution to modes derived at time  t  = 0 . 25  by ILDManalysis. The model parameters are  S  0 = 100 ,  E  0 = 1 ,  k  1 = 1 ,  k  − 1 = k  2 = 100 . with three dimensionless parameters:  =  E  0 S  0 ,  S  t  = k  2 k  − 1 ,  Q = K  m S  0 = k  − 1 + k  2 k  1 S  0 Our methods are not restricted to the local analysis near the steadystate, so we can investigate the dynamical interaction betweenspecies in the modes additionally to the analytical results of Palssonand Lightfoot (1984). 3.1.1 ILDM analysis in COPASI   Following the approach of Palsson and Lightfoot, we consider five possible combinationsof dimensionless parameters. In all considered cases, the cleartime-scale separation occurs.(i)   → 0 :  this is the classical situation of Michaelis–Mentenkinetics. In this limit, the motion on the fast time scale isdominated by the complex  C   decoupled from the substrate S  . On the slow part, the changes of substrate and complexare balanced (Palsson and Lightfoot, 1984). Here, the QSSassumption is valid.Our time-dependent classical ILDM leads to the distinctionbetween slow and fast modes (for chosen  tol = 10 − 3 ) attimes  t  > 0 . 001. The analysis of the transformation matrix T − 1 shows that the complex  C   dominates the fast mode.The contributions of both metabolites to the slow mode arecomparable and so, the changes of substrate and complex arebalanced [see Table 1, which displays the matrix ‘Modes( i ,  j )’at time point  t  = 0 . 25].The analysis via the modified method leads to theclassification of the complex  C   as a fast metabolite (however,for time  t  > 0 . 017 first). Thus, the QSSA for metabolites  C   is justified in this case. The last method does not indicate theinteraction between metabolites in the slow mode (compareresults in Table 1).(ii)  Q → 0(Fig.1):thislimitisveryinterestinginordertoconsiderthedynamicalinteractionbetweenmetabolitesinslowandfastmodes. For our analysis in COPASI, we choose   = 0 . 5 ,  S  t  = 1 ,  Q = 0 . 05 and tolerance is  tol = 10 − 2 .No reduction is possible by modified ILDM method until t  = 2 . 6. For  t  > 2 . 6, the substrate is distinguished as fast,its contribution to the fast space is  > 95%. The analysisof the classical ILDM shows that after the short transientphase ( t  < 0 . 025) the complex  C   dominates the fast space.This dominance decreases with time and for  t  ≈ 1 . 13, bothsubstrateandcomplexarebalancedalsointhefastmode(theircontributionstothefastspaceareequal).Thenthecontributionof the substrate increases and finally, the substrate dominatesthe fast space near the equilibrium. Contributions of both 2819   b  y g u e  s  t   on O c  t   o b  e r 1  6  ,2  0 1 4 h  t   t   p :  /   /   b i   oi  nf   or m a  t  i   c  s  . oxf   or  d  j   o ur n a l   s  . or  g /  D o wnl   o a  d  e  d f  r  om 
We Need Your Support
Thank you for visiting our website and your interest in our free products and services. We are nonprofit website to share and download documents. To the running of this website, we need your help to support us.

Thanks to everyone for your continued support.

No, Thanks