BISIMULATIONS FOR DIFFERENTIAL EQUATIONS - PowerPoint PPT Presentation

MIRCO TRIBASTONE IMT SCHOOL FOR ADVANCED STUDIES LUCCA JOINT WORK WITH LUCA CARDELLI, MAX TSCHAIKOWSKI, AND ANDREA VANDIN BISIMULATIONS FOR DIFFERENTIAL EQUATIONS http://sysma.imtlucca.it/tools/erode/

MOTIVATION 1 2 3 V in 1 2 3 Fig. 2: Mesh networks adapted from [26]. (Bio-)chemistry Systems biology Engineering ▸ Ordinary differential equations as a universal mathematical modelling language ▸ System complexity leads to large-scale models ▸ Reduction/abstraction needed to gain physical intelligibility and ease analysis cost


 
 
 
 POLYNOMIAL ODES ▸ Polynomial ODEs are the class of ODEs where the derivatives are multivariate polynomials (over the system’s variables) ▸ They cover mass-action kinetics , but they can also encode many nonlinearities (trigonometric, rational, exponential functions) ▸ For biochemistry polynomial ODEs may encode other kinetics such as Hill, Michaelis-Menten, sigmoids, etc. 
 x = 1 ˙ z := sin( x ) , w := cos( x ) x = 1 ˙ y = z ˙ y = sin( x ) ˙ z = w ˙ w = − z ˙

LUMPING DIFFERENTIAL EQUATIONS ▸ ODE lumping means finding a partition of variables such that each partition block (equivalence class) can be associated with a single equation ▸ The lumped ODE preserves the original dynamics: ▸ Forward lumping preserves sums of the solutions of the variables in each block ▸ Backward lumping identifies variables with the same solution in each block ▸ ODE lumping can be exact or approximate ▸ ODE lumping is complementary to other techniques, e.g. those for fast-slow decomposition (QE/QSSA) and domain-specific reductions such as fragmentation

MOTIVATIONAL EXAMPLE: MARKOV CHAIN LUMPING STRUCTURE DYNAMICS LUMPED DYNAMICS π 1 = − ( k 1 + k 2 ) π 1 ˙ 1 ˙ Π 1 = − ( k 1 + k 2 ) Π 1 k 1 k 2 π 2 = k 1 π 1 ˙ ˙ Π 2 = ( k 1 + k 2 ) Π 1 2 3 π 3 = k 2 π 1 ˙ q 2 , 1 = q 3 , 1 = 0 Π 2 ( t ) = π 2 ( t ) + π 3 ( t ) ▸ CTMC lumping as a special case of ODE lumping (for a class of linear ODEs) ▸ Can we generalise to nonlinear (polynomial) ODEs? ▸ What is the structural analogous to a transition matrix?

MASS-ACTION LUMPING AT A GLANCE x 1 = − x 1 + x 2 − 3 x 1 x 3 + 4 x 4 ˙ y 1 y 2 y 3 y 1 = − 3 y 1 y 2 + 4 y 3 ˙ x 2 = + x 1 − x 2 − 3 x 2 x 3 + 4 x 5 ˙ � { x 1 , x 2 } , { x 3 } , { x 4 , x 5 } y 2 = − 3 y 1 y 2 + 4 y 3 ˙ x 3 = − 3 x 1 x 3 + 4 x 4 − 3 x 2 x 3 + 4 x 5 ˙ y 3 = +3 y 1 y 2 − 4 y 3 ˙ x 4 = +3 x 1 x 3 − 4 x 4 ˙ x 5 = +3 x 2 x 3 − 4 x 5 ˙ 2 2 1.8 1.8 x 1 y 1 x 2 1.6 1.6 y 2 x 3 y 3 1.4 1.4 x 4 x 5 1.2 1.2 x 1 +x 2 x(t) y(t) x 4 +x 5 1 1 0.8 0.8 0.6 0.6 0.4 0.4 0.2 0.2 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 Time Time

FROM POLYNOMIAL ODES TO REACTION NETWORKS ▸ Main idea: take each monomial appearing in the derivative and transform it into an edge of a (labelled) bipartite multigraph: a reaction n x p i Y x k = . . . + α ˙ i + . . . i =1 Stoichiometric 
 From continuous to discrete “Reaction rate” coefficient Species REACTANTS PRODUCTS α p 1 X 1 + . . . + p n X n → p 1 X 1 + . . . + p n X n + X k − − No physical meaning, used only for reasoning on equivalences

FORWARD EQUIVALENCE ▸ Forward equivalence generalises CTMC forward lumpability ▸ Main intuition (borrowed from process algebra): projected to a species, other reactants are partners/communication channels (Multiset) partners of X 1 α X 1 + p 2 X 2 + p 3 X 3 → . . . − − Candidate 
 equivalent 
 (Multiset) partners of Y 1 species β Y 1 + q 2 Y 2 + q 3 Y 3 → . . . − − ▸ Two species are equivalent if they have equal aggregate rate toward any block, with any possible partners


 
 FORWARD EQUIVALENCE, FORMALLY Multiset of reactants Multiset of products α ρ i is the multiplicity of X i − − → π ρ FLUX NET STOICHIOMETRY FORWARD RATE P X j ∈ G φ ( X i + ρ , X j ) ✓ P ◆ i ρ i X φ ( ρ , X i ) := α ( π i − ρ i ) fr ( X i , ρ , G ) := , [ ρ ]! := [ X i + ρ ]! ρ 1 , . . . , ρ n α − → π all ρ ▸ A partition of species if a forward equivalence if, for any two blocks and any two species it holds that 
 H, H 0 X i , X j fr ( X i , ρ , H 0 ) = fr ( X j , ρ , H 0 ) for all multisets partners ρ ▸ Characterisation result , extending previous work 
 [CONCUR’15,TACAS’16] Maximal aggregation of polynomial dynamical systems Luca Cardelli a,b,1 , Mirco Tribastone c,1,2 , Max Tschaikowski c,1 , and Andrea Vandin c,1 a Microsoft Research, Cambridge CB1 2FB, United Kingdom; b Department of Computing, University of Oxford, Oxford OX1 3QD, United Kingdom; and c Scuola IMT Alti Studi Lucca, 55100 Lucca, Italy Edited by Moshe Y. Vardi, Rice University, Houston, TX, and approved July 28, 2017 (received for review February 16, 2017) Ordinary differential equations (ODEs) with polynomial derivatives variables in a single block. Furthermore, the freedom in choosing are a fundamental tool for understanding the dynamics of systems an arbitrary initial partition is instrumental to producing reduc- tions that preserve the dynamics of desired original variables, across many branches of science, but our ability to gain mecha- which are then not aggregated. nistic insight and effectively conduct numerical evaluations is crit- ically hindered when dealing with large models. Here we propose Mathematically, our approach is a generalization of well-


 
 
 BACKWARD EQUIVALENCE ▸ Analogous of backward lumpability in Markov chains: equivalent species have the same solution at all times ▸ Definition in the same bisimulation style, with a twist: 
 br ( X i , M , H ) = br ( X j , M , H ) where is an equivalence relation on multisets of species M naturally induced by the equivalence over species, e.g.: 
 A ∼ B = ⇒ A + B + C ∼ M 2 B + C ▸ Characterisation result: extension of [CONCUR’15] Maximal aggregation of polynomial dynamical systems Luca Cardelli a,b,1 , Mirco Tribastone c,1,2 , Max Tschaikowski c,1 , and Andrea Vandin c,1 a Microsoft Research, Cambridge CB1 2FB, United Kingdom; b Department of Computing, University of Oxford, Oxford OX1 3QD, United Kingdom; and c Scuola IMT Alti Studi Lucca, 55100 Lucca, Italy Edited by Moshe Y. Vardi, Rice University, Houston, TX, and approved July 28, 2017 (received for review February 16, 2017) Ordinary differential equations (ODEs) with polynomial derivatives variables in a single block. Furthermore, the freedom in choosing are a fundamental tool for understanding the dynamics of systems an arbitrary initial partition is instrumental to producing reduc- tions that preserve the dynamics of desired original variables, across many branches of science, but our ability to gain mecha- which are then not aggregated. nistic insight and effectively conduct numerical evaluations is crit- ically hindered when dealing with large models. Here we propose Mathematically, our approach is a generalization of well-

MAXIMAL AGGREGATION THROUGH PARTITION REFINEMENT FORWARD RATE BACKWARD RATE fr ( X i , ρ , H 0 ) = fr ( X j , ρ , H 0 ) br ( X i , M , H ) = br ( X j , M , H ) ▸ Both in the same style of probabilistic bisimulation (where the partners are analogues of action type) ▸ Intuition for computing the maximal aggregation through a partition refinement algorithm ▸ Polynomial time and space complexity (in the number of species, number of monomials and maximum degree) ▸ Extensions of the works of Derisavi et al. , Valmari & Franceschinis , Baier et al. , our own [MFCS’15,TACAS’16]

PARTITION REFINEMENT EXAMPLE Reaction Network ▸ Binding model k 1 A u,u → A p,u − k 2 A p,u → A u,u ▸ Occurs over one of two binding − k 1 sites when it is phosphorylated A u,u → A u,p − k 2 A u,p → A u,u − ▸ Classic, basic model in biochemistry k 3 A p,u + B → A p,u B − ▸ Intuition: if the binding sites are k 4 A p,u B → A p,u + B − identical then, by symmetry, explicit k 3 A u,p + B → A u,p B − identity is unimportant k 4 A u,p B → A u,p + B −

PARTITION REFINEMENT EXAMPLE B C Reaction Network Current partition k 1 � { A u,u , A p,u , A u,p , B, A p,u B, A u,p B } A u,u → A p,u − First iteration k 2 A p,u → A u,u − Set of splitters k 1 � A u,u → A u,p { A u,u , A p,u , A u,p , B, A p,u B, A u,p B } − k 2 sp A u,p → A u,u − fr ( A p,u , B, sp ) = − 3 k 3 A p,u + B → A p,u B − fr ( A u,p , B, sp ) = − 3 k 4 A p,u B → A p,u + B fr ( B, A p,u , sp ) = − 3 − k 3 fr ( B, A u,p , sp ) = − 3 A u,p + B → A u,p B − fr ( A p,u B, ∅ , sp ) = 4 k 4 A u,p B → A u,p + B − fr ( A u,p B, ∅ , sp ) = 4 ▸ The initial partition may be arbitrary ▸ Useful to single out observables that are not to be aggregated