A Model Checking Approach to Computational Systems Biology David Šafránek with Luboš Brim and Milan Češka Masaryk University Czech Republic FUNDP Namur, 9.5.2014 1/139 Q Introduction O LTL Model Checking 0 Discrete Abstraction of ODE Models Q Parameter Synthesis and Coloured Model Checking Q Case Studies • E. Coli Ammonium Transport • Cell Cycle Regulation Synthetic Biology: Trichlorpropane Degradation Q Parameter Synthesis and Classification for Boolean Networks FUNDP Namur, 9.5.2014 2/139 Introduction r____j^i___■ ____I r° _ I_______I iv a _ _i 9 E. Coli Ammonium Transport • Cell Cycle Regulation 9 Synthetic Biology: Trichlorpropane Degradation FUNDP Namur, 9.5.2014 Motivation Model Analysis in Systems Biology network reconstruction biological knowledge databases t model validation gene reporters, DNA microarray, mass spectrometry, ... Bacterial DNA biological network hypothesis emergent properties model questions model specification SBML, diferenciální rovnice, boolovská sít, Petriho sít, ... 1 model analysis static analysis, numerical simulation, analytical methods, model checking hypothesis testing, property detection, new hypothesis inference FUNDP Namur, 9.5.2014 4/139 Systems View of Processes Driving the Cell FUNDP Namur, 9.5.2014 5/139 Systems View of a Cell: Biological Networks • identify substances (macromolecules, ligands, proteins, genes, .. .) • identify interactions ((de)complexation, (de)phosphorylation, .. .) FUNDP Namur, 9.5.2014 6/139 Systemic View of the Cell: Biological Networks FUNDP Namur, 9.5.2014 7/139 Systemic View of the Cell: Biological Networks FUNDP Namur, 9.5.2014 7/139 Systems Biology of a Cell Literature System boundary metabolite Metabolic networks Signalling Pathways Protein-protein Interaction protein 1 A Metabolome space metabolite3 B rend a ENZYME Pr< teome space SwissProt ^_-J ranscriptome space Gene RNA1 Regulatory Pathways | ------- gene3 Genome space gene GenBank ■DDBJ. EMBL S^Ensemfih üernTme Browser FUNDP Namur, 9.5.2014 8/139 Biological Networks and Pathways 9 what is the "right" meaning? • in order to analyse we need to formalise FUNDP Namur, 9.5.2014 9/139 Biological Networks and Pathways Graphical Specification in SBGN Systems Biology Graphical Notation PD ER • PD: biochemical interaction level (the most concrete) • ER: relations among components and interactions • AF: abstraction to mutual interaction among activities FUNDP Namur, 9.5.2014 Graphical Specification in SBGN Systems Biology Graphical Notation • SBGN.org iniciativě (from 2008) 9 standard notation for biological processes • http://sbgn.org • Nature Biotechnology (doi:10.1038/nbt.l558, 08/2009) • three sub-languages: • SBGN PD (Process Description) (doi:10.1038/npre.2009.3721.1) • SBGN ER (Entity Relationship) (doi:10.1038/npre.2009.3719.1) • SBGN AF (Activity Flow) (doi:10.1038/npre.2009.3724.1) • tool support: • SBGN PD supported by CellDesigner • SBGN-ED (http://www.sbgn-ed.org) FUNDP Namur, 9.5.2014 12/139 Kinase Cascade in CellDesigner (SBGN) FUNDP Namur, 9.5.2014 Why to model? Thanateßh^JCte Dysplasia ■ short long bones • macrocephaly • low nasal bridge - spinal stenosis • temporal lobe malformations Nat Genet 1995, 9:321-8. e.g., FGFR3-related skeletal dysplasia FUNDP Namur, 9.5.2014 14/139 Why to model? Need to explain... FGF2 J FGFR3 Erk FGF2 lh» P. KrejCf, Masaryk University FUNDP Namur, 9.5.2014 Why to model? Knowledge is increasing... 2001 2012 FGFR3 STATI '/ CKI IL6, LIF, IL11, IFNv FGFR3 SOCS1/3 PKC STAT 1/3 STATI FrS2, Gab1, SHC proliferation CNP NPR-B cGMP RasŕRafíMEK/Erk P KG growth arrest growth arrest CKI MMP J matrix degradation FUNDP Namur, 9.5.2014 16/139 Model of a Biological System complex biological system g in vitro/in vivo f N formal model ^ = —ki ■ A + k2 ■ B M f = h-A-k2-B in silico -\-w- model parameters • model is an approximation of the biological system • built on first-principles and known hypotheses ^> e.g., elemental reactions, experimental observations, . .. • model is parametrized • parameters provide a space for refinement ^> typically quantitative information (e.g., reaction rate) FUNDP Namur, 9.5.2014 17/139 Systems View of the Cell 9 syntax of the systems model is a network: 9 components (nodes) - e.g. chemical substances • interactions (edges) - e.g. chemical reactions • each component is assigned some quantity • discrete: number of molecules • continuous: molecule concentration in a compartment (solution) • can be visualized by color intensity of a node • dynamics drives the change of node colour intensity in time • driven by global rules (e.g., mass-action reactions) FUNDP Namur, 9.5.2014 18/139 Biological Model Formal Definition Denote St = Z domain of stoichiometric coefficients. Biological model is a tuple (S, /?, reanet, regnet, map), where: • S C N ... (finite) species index set • R C N ... (finite) reactions index set • reanet C S x R ... reaction network • regnet CSxRx {inh, act} ... regulatory network o map : reanet —)► §t ... stoichiometric map Members of S are denoted: si, s^,.... Members of R are denoted: ri, r2,.... FUNDP Namur, 9.5.2014 19/139 Biological Model - Example FUNDP Namur, 9.5.2014 20/139 Biological Model - Example • S = {A1,A2,A3,A4,A5} « R = {n, r2, r3, r4, r5} • reanet, map: (ri) 24i + A2 -> /43 + /44 + A> (r2) /A4 + As -> 3/\2 (r3) -> /Ai (r4) A -> (r5) A3 -> FUNDP Namur, 9.5.2014 21/139 Biological Model - Example FUNDP Namur, 9.5.2014 Biological Model - Example • S = {A1,A2,A3,A47A5} 9 R = {ri, r2, r3, r4, r5} • reanet,map: (ri) 2AX + A2^ A3 + A4 + A5 (r2) /A4 + A5 3A2 (f3) -+ *i (r4) /Ai -)■ (r5) A3 ->• • regnet : A2 inhibits A3 activates fy, ^5 inhibits r\ FUNDP Namur, 9.5.2014 23/139 Model Semantics Discrete Case A + B->AB o state configuration captures number of molecules (#[AB],#[A],#[B])eN30 global rule: » one molecule AB is added to the solution • one molecule A is removed from the solution • one molecule B is removed from the solution #[AB](t + 1) = #[AB](t) + 1 #H(t + 1) = #[/\](t) - 1 #[B](t + i) = #[e](t)-i FUNDP Namur, 9.5.2014 Model Semantics Discrete Case Consider three reactions: A^ B A + B^AB AB^A + B 9 state configuration has the form (#/4, #B, #/46) e Nq o consider, e.g., configuration (2,2,1) =4> what is the next configuration? • reactions run in parallel ... FUNDP Namur, 9.5.2014 25/139 Model Semantics Continuous Case A + B^AB • continuous state captures concentration of molecules in a certain volume (the solution): ([AB},[Al[B])eR3+ global rule: • a mass of AB outflows from the solution • a mass of A inflows to the solution • a mass of B inflows to the solution d[AB] dt = V d[A] _ d[B] _ dt ~ dt ~ v where v = /c[/4][B], k is the reaction rate constant. The law of mass action. FUNDP Namur, 9.5.2014 26/139 Model Semantics Discrete Gene Regulatory Networks A e {0,1,2}, ße{0,l} tAA = 2, tAB = 1 Ka,9 = 2 Ka,{a} = 0 Kr fll = 0 K e,0 B,{A} = 1 0,0 I 0,1 T 2,1 • introduced by René Thomas [1973] • refined by Chaouiya et al. [2003] FUNDP Namur, 9.5.2014 Model Semantics species S are interpreted as model variables • boolean models: val(Si) £ {present, absent} • discrete-value models: val(S;) E No • continuous-value models: val(Si) £ current values of all model variables make the state reaction is interpreted as a rule that affects (changes) the state Variables are always considered bounded (maximal values can be given by physical limits, e.g., the cell volume). FUNDP Namur, 9.5.2014 28/139 Rule Interpretation Modelling of time • exact time of reaction occurence =4> continuous-time models • time of reaction occurence abstracted =4> discrete-time models (ticked or untimed) Modelling of noise • deterministic rules - noise absent (large populations) =4> always one possible execution under the same conditions • stochastic rules - noise present (small populations) =4> many different executions possible under the same conditions FUNDP Namur, 9.5.2014 29/139 Model Semantics Spectrum - Brief ■ quantitative parameters ignored quantitative parameters required I ] variables -s=- -3*- discrete continuous I FUNDP Namur, 9.5.2014 30/139 Model Semantics Spectrum - Detailed & ^ ^ * co CD > time qualitatively abstracted ^ CD -I—« CD l_ O CO =3 time quantitatively modeled v& o V "■4—' o o XT FUNDP Namur, 9.5.2014 31/139 From Biological Hypotheses to Temporal Properties a wet-lab measurements =4> time-series data 9 low resolution - e.g., microarray data, series of western blots • high resolution - fluorescence measurements (e.g., gene reporters) • most typically population-level measurements (average behaviour) • literature provides other constraints on system dynamics • e.g., multiple steady states, species concentration correlation, • all can be formally captured by means of temporal logics FUNDP Namur, 9.5.2014 32/139 Wet-lab Measurements 9 systems measurements of transcriptome (mRNA concentration) • very imprecise! FUNDP Namur, 9.5.2014 33/139 Wet-lab Measurements • western blots • measurements of protein binding (presence of certain proteins) FGF2 CI y ]»' MV Ih 2h 41i C2 lP:Frs2 t~ Shp2 — Grb2 Frs2 _ IP; Gabl ~Shp2 WB Gib2 ---- Gab) \p Slit Shp2 ■ WD r^6 She P52 FUNDP Namur, 9.5.2014 34/139 O Model is built on first-principles =4> purely qualitative (network topology) =4> quantitative aspects represented by parameters O To build an executable model we need to find all possible constraints that can be formulated. =4> static and dynamic constraints (properties) O To find admissible parameter values we need further elaboration =4> fitting to wet-lab measurements is a problem when some data are too imprecise (practical identifiability) FUNDP Namur, 9.5.2014 35/139 Qualitative vs. quantitative temporal properties • qualitative properties (LTL, CTL) • modalities (possibilities/necessities in future behaviour) • reachability of particular (sets of) states • temporal ordering of events, monotonicity • temporal correlations of model variables • stability (attractors, basins of attraction) 9 quantitative properties • deterministic (MTL, MITL, STL) 9 enhance modalities with (dense) time information o exact timing of events, time-bounds o stochastic (PLTL, PCTL, CSL) • probability of property satisfaction • stochasticity combined with time FUNDP Namur, 9.5.2014 36/139 Qualitative vs. quantitative temporal properties • qualitative properties (LTL, CTL) • modalities (possibilities/necessities in future behaviour) • reachability of particular (sets of) states • temporal ordering of events, monotonicity - time-series • temporal correlations of model variables - time-series • stability (attractors, basins of attraction) 9 quantitative properties • deterministic (MTL, MITL, STL) 9 enhance modalities with (dense) time information o exact timing of events, time-bounds o stochastic (PLTL, PCTL, CSL) • probability of property satisfaction • stochasticity combined with time FUNDP Namur, 9.5.2014 36/139 Temporal Property Examples Qualitative properties enzyme E is never permanently exhausted GF(E > 0) • all molecules of the substrate S are finally transfered to the product P provided that the final state is stable S == 5 FG(P == 5 A S == 0) • enzyme E is used and finally returned back (E > 2) U [(0 < E < 2) U (E > 2)] FUNDP Namur, 9.5.2014 37/139 Temporal Property Examples Quantitative properties • in the first 10 time units, enzyme E cannot permanently exhausted G[°'10lF(E>0) q all molecules of the substrate S are transfered to the product P minimally in 2 and maximally in 5 time units provided that the final state is stable S == 5 FP'5]G(P == 5 A S == 0) • enzyme E is used and finally returned back within the given time intervals (E > 2) Ul1'2' [(0 < E < 2) IP'2! (E > 2)] FUNDP Namur, 9.5.2014 38/139 Temporal Property Examples • oscillation LTL: (G[(A < 3) F(A > 3)]) A (G[(A > 3) F(A < 3)]) • bistability CTL: EFAG(4 < 5) A EFAG(4 > 5) • probabilistic modality PCTL: P>0.g[F(>4 = 3)] • probabilistic modality with time CSL: P>o.9[F[1-2](^ = 3)] FUNDP Namur, 9.5.2014 39/139 System Construction and Formal Methods required properties verification specified properties FUNDP Namur, 9.5.2014 40/139 Knowledge Discovery and Formal Methods hypothesis validation prediction inferred properties FUNDP Namur, 9.5.2014 I I I kl VUU^tlVI I O LTL Model Checking %SV LylbLictc MUbT-i dLT-IUll UT v>/|_>/C IVIUUclb • E. Coli Ammonium Transport • Cell Cycle Regulation 9 Synthetic Biology: Trichlorpropane Degradation FUNDP Namur, 9.5.2014 Definition Let AP be the set of atomic propositions (logical expressions over model variables, typical inequalities). Kripke structure is the quadruple K = (S, So, T, L) where: • S is the finite set of states • So Q S is the set of inititial states • T C S x S such that Vs e S, 3s; e S : (s, s;) e 7" • L is the labeling L : S —>► 2^p FUNDP Namur, 9.5.2014 Kripke structure - properties • for a state s e S, L(s) represents the set of all atomic propositions satisfied in s • unfolding of the Kripke structure from any initial state is always an infinite-depth tree 9 maximal paths in the unfolding represent individual (infinite) executions of the Kripke structure FUNDP Namur, 9.5.2014 44/139 Linear-time Temporal Logic - syntax Let AP be the set of atomic propositions. Formula cp is linear temporal logic (LTL) formula iff the following holds: • cp = p for any p 6 AP 9 If cpi and LTL formulae then: • ""^l, <£i A y?2 and y?i V p2 are LTL formulae • Xpi, F(pi a Gy?i are LTL formulae • piUp2 a piRp2 are LTL formulae FUNDP Namur, 9.5.2014 45/139 Linear Temporal Logic - semantics Let 7r = so, si,s/,... be an infinite sequence of states (a path) in a Kripke structure K. For j > 0 we denote 7T7 the suffix s/, sy+i,s/,.... Satisfiability relation |= is defined by induction: • 7r |= p iff p e Z.(sb) • 7T |= \ff 7T ty= if 9 7T \= (fl A 2 • 7T |= X(fi iff 7T1 |= o 7r |= F(p iff 3/ > 0.7r' |= (fi O 7T |= Gy? iff V/ > 0. tv1 |= • 7T |= ^iU(^2 iff 3/ > 0- ^ |= ^2 and V/ < J. tv1 |= <£i 9 7T |= (fiR(f2 iff V/ > 0, V0 < / < j. 7T1 ^ (fil ^ 7TJ \= (f2- FUNDP Namur, 9.5.2014 46/139 Linear Temporal Logic - semantics FUNDP Namur, 9.5.2014 47/139 Linear Temporal Logic - semantics For any formulae cpi, cp2 the following holds: —if cp = G—ap The full expressiveness is achieved by using just the operators n,A,X,U. FUNDP Namur, 9.5.2014 48/139 Linear Temporal Logic - semantics For any formulae cpi, cf2 the following holds: —iFcp = G^(p -i(^lU^2) = -«^lR-"^2 The full expressiveness is achieved by using just the operators -,A,X,U LTL formulae are most typically interpreted universally over Kripke structure paths: FUNDP Namur, 9.5.2014 48/139 Linear Temporal Logic - semantics For any formulae tpi, tp2 the following holds: -iFy> = G—ap -"(^iU^) = —«piR—«p2 The full expressiveness is achieved by using just the operators -,A,X,U. LTL formulae are most typically interpreted universally over Kripke structure paths: Kripke structure as a model for a formula Let K be a Kripke structure. A formula tp is satisfied by K, K |= cp iff for each execution tt = Sq, ... such that Sq G Sq it holds tt |= A, the set of all (infinite) words for which there exist a corresponding accepting run of A, 9 cj-regular languages are closed under complementation. FUNDP Namur, 9.5.2014 52/139 Biichi automata examples true FUNDP Namur, 9.5.2014 53/139 LTL Model Checking 9 Specification formalized as LTL formula Automata-based approach to LTL model checking 9 Employs Buchi automata to express • all paths of the Kripke structure under consideration • all paths violating the specification • Model satisfies the specification if the intersection of the sets is empty, i.e., if the synchronized Buchi automaton accepts empty language. LTL model checking problem is reduced to the detection of accepting cycles in the graph of a Buchi automaton. FUNDP Namur, 9.5.2014 Model Checking as a language inclusion problem Interpretation of a path tv = s&,si,... in a Kripke structure K is a sequence of sets of APs: L(tt) = /.(sö),/.(si),... Problem For a given Kripke structure K = (S, So, T, Z.) and a given LTL formula p decide K |= p. Reformulation Let Z = 2AP. Consider two languages of infinite words: O Ck — {L{tt) I 7r is a path in K} O = {L(7r) I 7T |= V9} Then K \= cp iff £k C FUNDP Namur, 9.5.2014 55/139 Kripke structure as a Büchi automaton Claim For each Kripke structure K = (S, So, 7~, L) we can construct a Buchi automaton Ak such that Ck = C{A^)\ • AK = (S,2AP,S0,6,S) where q e 5(p, a) 44> (p, q) e 7" A L(p) = a. Observation Note that F — S (the set of final states coincides with the state space). FUNDP Namur, 9.5.2014 56/139 LTL formula as a Biichi automaton Theorem [Vardi, Wolper 1986] For each LTL formula cp there exists (and can be effectively constructed) a Biichi automaton such that = £(A^)- Construction goes through a generalized BA (extended in the acceptance condition - a system of accepting states sets, requirement to infinitely often visit all accepting sets). Complexity is 2°(n^ where n is the size of the formula. There exist many algorithms - check, e.g., http://spot.lip6.fr/wiki/. Note LTL is less expressive then BAs. FUNDP Namur, 9.5.2014 Let A = (Sa,T,S0a,Sa,Sa), B = (Sb, S0e, Ffi) be Buchi automata, and FA — SA. Then a Buchi automaton Ax B that accepts the language L(A x B) — L(A) n /-(B) can be constructed in the following way: • A x B = (SA x Se,5I,So4 x S0e,^x8,S^ x F8), • (p', q;) G 5AxB((p, q), a) for all p' e 5A(p, a) and q' e 5B(q,a). FUNDP Namur, 9.5.2014 Model Checking reduced to language emptyness problem Claim For each LTL formula cp it holds that co-C^A^) = C(A^) • K • K • K • K • K p Ck Q £cp P & L(AK) C L{AV) ip & L(AK) n co-C(Ap) = 0 p L{AK) n L(A^) = 0

AA V/ G {1,n}. V/ G {1,m}. crt c/X/ FUNDP Namur, 9.5.2014 67/139 Reaction Kinetics • format of chemical reactions: 7i*i + '' -+7mXm -^S1Y1 + " - + 5nYm ri e {0,e N note we expect {Xi,Xm} n { Yi,Yn} = 0 o subclass of general mass action kinetics: g(Xll...,Xm) = 5ikX^X^---XZm g(Xi,Xm) = -r.kXfX*2 ■■■xi- • corresponds to the class of multi-affine autonomous systems • limitation: homodimerization A + A —>> AA • reactions of the form X SiYi + - — + 5nYm5j e N result in affine systems V/ e {1,r?}. V/ G {1,rn}. jv2 eft C/X; FUNDP Namur, 9.5.2014 67/139 Regulatory Kinetics protein dynamics driven by protein-regulated transcription Hill kinetics approximated in terms of ramp functions X —Y dY ~dt = Ar+(X,x£i,x£2) 9 k G M+ is kinetic parameter r+(X,xtl,xt2) rxi ramp functions can describe cooperative regulations by means of summation and multiplication FUNDP Namur, 9.5.2014 68/139 General Kinetic Models • both kinetics combined 9 right-hand side of any ODE is a mapping g"(x,p) where p is a vector of unknown parameters • (piece-wise) multi-affine in x • affine in p o these properties enable us to (are necessary to): • make a discrete finite overapproximation of the system dynamics • discretize the parameter space - possible values of p synthesis of unknown parameters FUNDP Namur, 9.5.2014 69/139 Rectangular Abstraction for Kinetic Models Overapproximative Abstraction on Rectangles FUNDP Namur, 9.5.2014 70/139 Rectangular Abstraction for Kinetic Models approach of [Belta, Habets, van Schuppen] continuous phase-space is bounded and abstracted by non-deterministic automaton B • r- A FUNDP Namur, 9.5.2014 Rectangular Abstraction for Kinetic Models approach of [Belta, Habets, van Schuppen] continuous phase-space is bounded and abstracted by non-deterministic automaton K K—K B A FUNDP Namur, 9.5.2014 Rectangular Abstraction for Kinetic Models approach of [Belta, Habets, van Schuppen] continuous phase-space is bounded and abstracted by non-deterministic automaton K—K B I- I A FUNDP Namur, 9.5.2014 Rectangular Abstraction for Kinetic Models Partition Definition Let X C Mn be a closed full-dimensional polytope. A partitioning Xpart(X) = {X/ I / = 1,..., m} of X is called admissible if O for all / = 1,..., m: X; is a closed full-dimensional polytope in Rn, e u^X; = x, O for all ij = 1,..., m, i ^ j, the intersection X/ n Xj is either empty, or a common face of X/ and Xj. If both polytope X and all subpolytopes X/, (/ = 1,..., m) are r?-dimensional rectangles, then an admissible partitioning Xpart(X) is called rectangular. FUNDP Namur, 9.5.2014 72/139 Rectangular Abstraction for Kinetic Models Piecewise Affine and Multi-Affine Mapping Definition A mapping g : X —> Rn is called piecewise-affine on Xpart(X) if the following two conditions hold: O g is continuous on X O for all / = 1,..., m there exist A\ G Rnxn and a; G Mn such that for all x G X/: g"(x) = /A/x + a,-, i.e. g- |x;- is an affine mapping. A mapping g : X —>» Rn is called multi-affine on Xpart(X) if the following two conditions hold: Q g is continuous on X O for all / = 1,..., m, g |x, is multi-affine, i.e. g |x, is affine w.r.t. every of its variables, while keeping all other variables constant. FUNDP Namur, 9.5.2014 73/139 Rectangular Abstraction for Kinetic Models Definition A piecewise-affine system on a polytope is a tuple X = (X,Xpart(X),Xb,to,g), where state set X is a full-dimensional polytope in IR", Xpart(X) is an admissible partitioning of X, xq e X is the initial continuous state, to £ R is the initial time, and g : X —>> IRn is a piecewise-affine function on Xpart(X). A trajectory x : [to, ti] —>• X of system x is a solution of the differential equation x(t) = g(x(t)), x(t0) =x0, (1) where ti is either the time instant that the trajectory leaves state polytope X, or t\ — oc, if trajectory x(t) remains in X forever. FUNDP Namur, 9.5.2014 74/139 Consider the affine system Z on rectangle [0,2] x [0,2] given by Obviously, (1.7,1.3)T is the unique steady state of this system. We partition the state set X into four squares: X(o,o) = [0,1] x [0,1], X(1)0) = [1,2] x [0,1], X(0)i) = [0,1] x [1,2], X(M) = [1,2] x [1,2]. FUNDP Namur, 9.5.2014 75/139 Rectangular Abstraction for Kinetic Models One may distinguish systems with the same dynamics on all polytopes in the partitioning, and systems with different dynamics on each subpolytope. In the second case, the dynamics on the boundary of two polytopes is still assumed to be continuous. FUNDP Namur, 9.5.2014 76/139 Rectangular Abstraction for Kinetic Models One may distinguish systems with the same dynamics on all polytopes in the partitioning, and systems with different dynamics on each subpolytope. In the second case, the dynamics on the boundary of two polytopes is still assumed to be continuous. Definition If X is an r?-dimensional rectangle, Xpart(X) is a rectangular partioning of X, and g : X —>> IRn is multi-affine on Xpart(X), then X = (X, Xpart(X), xq, to,g) is called a multi-affine system on rectangles. FUNDP Namur, 9.5.2014 76/139 Rectangular Abstraction for Kinetic Models Exit Facets Exit Facet Let x — 0^? XpaTt(X),xo, to,g) be a piecewise-affine system on a polytope X. A facet F of subpolytope X; is called an exit facet if there exists a trajectory of system Z, starting in X/, that attempts to leave X/ in finite time by crossing facet F. FUNDP Namur, 9.5.2014 77/139 Rectangular Abstraction for Kinetic Models Exit Facets Exit Facet Let x — (X, Xpart(X), xn, to,g) be a piecewise-affine system on a polytope X. A facet F of subpolytope X/ is called an exit facet if there exists a trajectory of system Z, starting in X/, that attempts to leave X/ in finite time by crossing facet F. Observation Let np denote the normal vector of F, pointing out of subpolytope X/, and let the affine dynamics on X/ be described by x = A\x-\- a\. Then F is an exit facet if and only if there exists x 6 F such that ^(Ax + a/) > 0. (2) Since the dynamics x = A\x + a; is affine, it suffices to check condition (2) on V(F), i.e. on the set of all vertices of facet F. FUNDP Namur, 9.5.2014 Rectangular Abstraction for Kinetic Models Exiting a polytope On which facet the trajectory exits a polytope X;? • if the trajectory leaves X/ through a point in the relative interior of a facet F, then it continues to an adjacent polytope Xj such that X/ n Xj = F, • if it leaves through a point on a lower-dimensional face, a problem arises since the face can be shared by more than two polytopes =4> this possibility is excluded and considered as singular (it is replaced by a sequence of several adjacent transitions executed in the same time instant) FUNDP Namur, 9.5.2014 78/139 Rectangular Abstraction for Kinetic Models Exiting a polytope Problem On which facet the trajectory exits a polytope X;? J • if the trajectory leaves X/ through a point in the relative interior of a facet F, then it continues to an adjacent polytope Xj such that X/ n Xj = F, 9 if it leaves through a point on a lower-dimensional face, a problem arises since the face can be shared by more than two polytopes this possibility is excluded and considered as singular (it is replaced by a sequence of several adjacent transitions executed in the same time instant) The rectangular abstraction abstracts from time. ] FUNDP Namur, 9.5.2014 78/139 1 íl \ x \ A(o.i> \ \ ^ \ / / / \ / \ 0 1 2 xj FUNDP Namur, 9.5.2014 79/139 Rectangular Abstraction for Kinetic Models Exiting a polytope Problem Does the trajectory leave a polytope X; in finite time? Theorem [Habets, Collins, Schuppen 2006] Consider an affine system x(t) = A;x(t) + a; on a closed full-dimensional subpolytope X; C There exists an initial state xq e Xj such that for all times t e T = [to, oo) the state trajectory belongs to the subpolytope, i.e. x(t; to,xo) C X/ if and only if there exists a fixed state in subpolytope X/. FUNDP Namur, 9.5.2014 80/139 Rectangular Abstraction for Kinetic Models Exiting a polytope Problem Does the trajectory leave a polytope X; in finite time? Lemma Consider an affine system x(t) = A\x{t) + a; on a closed full-dimensional subpolytope X/ C There exists an x 6 X/ such that A;x + a; = 0 if and only if 0 e ConvexHull({Au + 3; \ v e V(X/)}), (3) i.e. if and only if the zero vector is a convex combination of the direction vectors at the vertices. • alternatively numerical approaches can be used FUNDP Namur, 9.5.2014 81/139 Rectangular Abstraction for Kinetic Models Exiting a polytope O Subpolytope X; contains a fixed point, and at all vertices of Xj, the direction vector of the differential equation is pointing inward. In this case all trajectories that enter subpolytope X; will remain in X/ forever. Q Subpolytope X; does not contain a fixed point. Then all trajectories that enter X/ leave X/ in finite time. O Subpolytope X/ contains a fixed point, and there exists a vertex of X; where the direction vector of the differential equation is pointing out of X/. I.e., there exist trajectories that leave X/ and also trajectories that do not FUNDP Namur, 9.5.2014 82/139 Rectangular Abstraction for Kinetic Models The Abstraction Let x — 0^? ^part(^), *o, to,g) a piecewise-affine system, N — Impart(^)|- We construct a Kripke structure Kx = (S, So, T, L) representing the rectangular abstraction of x- 9 S = {si,s/\i} and we define a bijective map n : Xpart(X) S such that l~l(X/) = s/, • So = {s;} such that xq G n_1(s/) and x(t; £o?xo) C n_1(s/) for all t G (to, to + e) for some e > 0 • (s/, s/) e T if there exists x G n_1(s/) such that g-(x) = 0 • for every facet F — X; n X/ for that there exists a vertex v G V(F) satisfying njg(v) > 0, (n(X/), n(X;)) G 7 FUNDP Namur, 9.5.2014 83/139 Rectangular Abstraction for Kinetic Models Extension to multi-affine systems Rectangular abstraction can be employed also for (piecewise) multi-affine systems (proved only for rectangular polytopes). C. Belta, L.C.G.J.M. Habets, and V. Kumar. "Control of multi-affine systems on rectangles with applications to hybrid biomolecular networks." In Proc. 41th IEEE Conf. on Decision and Control, pages 534-539, New York, 2002. IEEE Press. Problem A sufficient and necessary condition for exiting a rectangle in finite time is not known. Theorem Let x(t) = g"(x(t)) be a multi-affine system on an r?-dimensional rectangle /?/ C If there exists a vector i^/GR" such that for all vertices v 6 V(/?/) we have wTg(v) > 0, then all state trajectories of this system leave rectangle /?/ in finite time. FUNDP Namur, 9.5.2014 84/139 Rectangular Abstraction for Kinetic Models Let x — {X, Xpart(X), xq, to,g) a piecewise-affine (or piecewise multi-affine) system and Kx its rectangular abstraction. Global necessity If for every path tt = s&,... in Kx, so G So, there exists an initial point xo G n(s/) such that the trajectory x(t; to?xo) of X corresponds to ir, i.e., x C (Js.G7r(n_1(s/)). Global sufficiency If for every trajectory x = x(t; to5xo) of x there exists a path 7r = so,... for some so G So such that xq e l~l(so) and x corresponds to 7r. FUNDP Namur, 9.5.2014 Temporal Properties for the Abstraction Kripke Structure • reachability 9 global: regardless the initial state, B eventually falls below 2 • local: if B initally below 2 then B does not exceed 2 • temporal properties • there is no initial state from which A falls below 6 before A exceeds 6 10 6--\ ^—• "5 0 time progress • properties defined by cj-regular languages • many useful properties can be formulated in LTL • some properties may require branching time (e.g., reachability of multiple steady state) FUNDP Namur, 9.5.2014 86/139 Rectangular Abstraction and Model Checking Let Kx be a rectangular automaton for a system x that is either (piecewise) affine or (piecewise) multi-affine. Let cp be an cj-regular property. 9 global sufficiency holds 9 Kx |= (f =^> x preserves ip • global necessity does not hold • Kx ^ ip does not necessarily imply "x does not preserve ip" • there might exist paths in Kx for which there is no trajectory in S, the reasons are of two kinds: O the abstraction combines behaviour of different trajectories =^> in piecewise-affine and multi-affine systems O known condition for exiting a rectangle in finite time is not sufficient =^> in multi-affine systems P. Collins, L. Habets, J.H. van Schuppen, I. Cerna, J. Fabrikova, and D. Šafránek. "Abstraction of Biochemical Reaction Systems on Polytopes", In Proceedings of the 18th IFAC World Congress. IFAC, 2011. pages 14869-14875 FUNDP Namur, 9.5.2014 87/139 Other Approac • regulatory kinetics abstracted by step functions • results in a piecewise-affine abstraction with different dynamics on individual rectalngles • gives a qualitative abstraction that is an overapproximation of original system • faces must be also included in the abstraction, trajectories are not continuous on faces =4> large state spaces, more expensive successor function good representation of regulatory logic (the extent of overapproximation is reasonable) H. de Jong, J.-L. Gouze, C. Hernandez, M. Page, T. Sari, J. Geiselmann (2004), Qualitative simulation of genetic regulatory networks using piecewise-linear models, Bulletin of Mathematical Biology, 66(2):301-340. FUNDP Namur, 9.5.2014 88/139 From Non-Linear to Piecewise (Multi)-Affine • a large class of molecular mechanisms modeled at activity-flow level (e.g., signalling pathways, gene regulatory circuits, ...) o optimal approximation of sigmoid functions by piece-wise affine functions (ramps) [Grosu et al. CAV 2011] model concentration ot [A] abstraction kinetics piece-wise multi-affine transient over-approximated sigmoidal kinetics steady state over-approximated mass action piece-wise affine transient over-approximated first-order sigmoidal steady state exact kinetics FUNDP Namur, 9.5.2014 89/139 %SV LylbLictc MUbT-i dLT-IUll UT v>/|_>/C IVIUUclb Q Parameter Synthesis and Coloured Model Checking • E. Coli Ammonium Transport • Cell Cycle Regulation 9 Synthetic Biology: Trichlorpropane Degradation FUNDP Namur, 9.5.2014 Motivation: Dynamical Systems with Parameters influencing reality parameters biological reality CLr eco n st r u et i o rT> mathematical model C^Tparameter tunmgj^> CTo bse rva t i o rT^> I parameter synthesis observed properties formalisation specified properties required properties admissible parameters FUNDP Namur, 9.5.2014 91/139 Problem Formulation Parameter Synthesis behavior constraints Si--.. C ^ Model M{p) : * = f(x, P) parameter constraints 4>i Pi 0 Parameter Synthesis Problem Assume P is the admissible parameter space. Given a behaviour constraint p, parameter constraint 0/, and a parameterised model A4, find the maximal set P CV of parameterisations such that p |= 0/ and M(p) |= p for all pG P. behaviour constraints V J r A parameter constraints dTormal 1 isatiqru^ valid parameter valuations ODE model C^app roxi m a t i orT^> ♦ PWMA model C^abstra^tionJ^ parameterised Kripke structure FUNDP Namur, 9.5.2014 93/139 Work Chronology Related Work • Batt et al. 2007: RoverGene, BDD/Polytopes-based approach • Batt et al. 2010: GNA, symbolic approach, piecewise affine • Grosu et al. 2011: RoverGene revisited, approximation improved • Bogomolov et al. 2015, SpaceEx, multi-affine hybrid automata Our Contribution • HIBI 2010, TCCB 2012: coloured LTL model checking, piecewise multi-affine, parallel algorithm o CMSB 2015: coloured CTL model checking, piecewise multi-affine, parallel algorithm • parameters represented as intervals • limitation: independent parameters only • ATVA 2016, CMSB 2016: parameters represented in first order logic, SMT solver employed, interdependent parameters • HSB 2015, FM 2016: discrete bifurcation analysis by coloured CTL model checking FUNDP Namur, 9.5.2014 94/139 Step 3: Parameter Synthesis Phase Space Discretisation Leads to Parameter Space Discretisation B 5 dA Hi = -ki ■ A + k2- B dB ~di = ki- A - k2- B k2 = 0.8 /ci = 0.6 o 2.5 5 A FUNDP Namur, 9.5.2014 95/139 Step 3: Parameter Synthesis Phase Space Discretisation Leads to Parameter Space Discretisation B 5 dA Hi dB ~di = -ki ■ A + k2- B = ki- A - k2- B k2 = 0.8 /ci = 0.6 2.5 \\\\\ \ \ \ \ ^ \ \ \ \ ^ \ \ \ ^ * \ \ \ ^ ^ \ \ \ * * ■ sj * a J J j 4 -i v r ' S, \ Si ^ * ^ \ a a J ' i( ^ A A 1 j a -i r r r i-i r r * * K K K r r K N * - r K ^ \ N ■ K \ \ \ \ ^ \ \ \ \ ■s. \ \ \ \ 0 2.5 5 A FUNDP Namur, 9.5.2014 95/139 Step 3: Parameter Synthesis Phase Space Discretisation Leads to Parameter Space Discretisation B 5 dA Hi = -ki ■ A + k2- B dB ~di = ki- A - k2- B k2 = 0.8 /ci = 0.6 5 A FUNDP Namur, 9.5.2014 95/139 Step 3: Parameter Synthesis Phase Space Discretisation Leads to Parameter Space Discretisation B 5 dA Hi = -kx ■ A + k2- B dB ~di = ki- A - k2- B k2 = 0.8 /ci = 0.6 5 A FUNDP Namur, 9.5.2014 95/139 Step 3: Parameter Synthesis Phase Space Discretisation Leads to Parameter Space Discretisation (0,0.4) (0.4,0.8) (0.8,1.6) (1.6,max) o \ \ \ \ \ \ \ €> \ \ \ o \ \ \ 0 \ \ \ FUNDP Namur, 9.5.2014 95/139 Step 3: Parameter Synthesis Phase Space Discretisation Leads to Parameter Space Discretisation B 5 dA Hi dB ~di = -ki ■ A + k2- B = ki- A - k2- B k2 = 0.8 ki = ? 0 \ N \ \ \ \ \ \ \ ^ \ \ \\|; 5 A (0,0.4) (0.4,0.8) (0.8,1.6) (1.6,max) o \ \ \ \ \ \ \ €> \ \ \ o \ \ \ 0 \ \ \ -2.5 • fc > 0 -2.5-/ci+2.5-/c2 > 0 FUNDP Namur, 9.5.2014 95/139 Step 3: Parameter Synthesis Phase Space Discretisation Leads to Parameter Space Discretisation dA Hi = -kx ■ A + k2- B dB ~di = ki- A - k2- B k2 = 0.8 ki = ? B 5 5—S—s—3I—O I \ * " J :0-H->0 0 \ \ \ \ \ \ \ 5 A 0 stateOO—estate 10 = -2.5 • /ci > 0 V -2.5 • /ci + 2.5 • k2 > 0 The transition exists if and only if the formula is satisfiable. Local parameter constraints are predicates over reals. FUNDP Namur, 9.5.2014 95/139 Parameter Synthesis by Coloured Model Checking parameterized Kripke structure of the model [B] 2.5 £2. f-^1 ► ► 2.5 [A] o o CD Q. O < identify states and colors for which the property does/doesn't hold parameter intervals where the specification is guaranteed (some might be missing) parameter intervals where the specification might be violated FUNDP Namur, 9.5.2014 96/139 Parameterised Kripke Structures State Transition Systems with Parameters Transitions with Parameters (coloured transitions) o each parameter valuation represents one Kripke structure • shared state space, different transition space FUNDP Namur, 9.5.2014 97/139 Parameterised Kripke Structures State Transition Systems with Parameters Transitions with Parameters (coloured transitions) o each parameter valuation represents one Kripke structure • shared state space, different transition space FUNDP Namur, 9.5.2014 97/139 Model Checking of Parameterized Kripke Structures Idea • for a model M and finite parameter space P consider Km — {Pi $i 5o, T, L) a parametrized Kripke structure • represent each parameterization by a distinct colour pE? • assume all transitions for each parameterization adequately coloured • find accepting cycles and get colours enabling accepting runs Procedure O construct the parametrized product BA of Km and the property BA O compute initial mapping of colours to states (state coloring) ^> propagate colours through the entire graph (BFS reachability) states on accepting cycles know all colours by which they are reached Q for each reachable accepting cycle aggregate (scan) the valid colours FUNDP Namur, 9.5.2014 State Coloring Let V denotes the set of all parameterizations. Further let /C = (P, S,PxI, So, 5, F) a parameterized product BA and let ttjG S, P CV. Succ^y, P)(a) = {p G P | 7A+ a} VS' C S. Succ(S', P) = |J Succ(7, P) Initial coloring: Sl/cc(So,P) Transition-enabling colours: a4/5 denotes /3 G (p, /-(a))) where p G P, is omitted to simplify the notation. FUNDP Namur, 9.5.2014 99/139 State Coloring Computation Compute Succ(S', P) over the PKS fC: Require: K = {V,S,V x E,S0,5,F), PCV, S' CS Ensure: R[a] = Succ(S', P)(ot) l: for all a G S do 2: R[a] P,ae S'} 5: while C? ^ 0 do 6: remove (a, P) from Q 7: if P g then 8: R[a] f- /?[«] U P 9: (?^(?©{(/3,PnP(o,/3))|o^/3,/3eS} 10: end if 11: end while • Q(a) = {p G V | 3P C P. p G P A (a, P) G Q} • (?©(?' = {(a, P) | P = Q(a) U • get S E ai P.. nd n .CO 7> 1 lours is reduced (levels of BFS) Chal lenges l number of states exponential w.r.t. number of variables • size of the parameter space exponential w.r.t. number of unknown parameters • many computations performed on a single graph FUNDP Namur, 9.5.2014 102/139 Parallel Implementation o multi-core data-parallel implementation of colour mapping propagation 9 states evenly distributed among threads by a hash-function 9 each thread responsible for a unique partition of colour mapping o threads communicate via a colour mapping update qeue (Q) o implemented as a set of lock-free qeues • one qeue per thread • threads synchronize on BFS levels FUNDP Namur, 9.5.2014 103/139 I I I I- I VJ %mA L> I V-/ I_ I I_ IVI \JUCI v_» I ICLr\11 I^ J^E^^ I ^ I c v es ~\~ es /\ ^ I f \ -f I ^ ^_ |\ yi ^ I ^ w£m UlbLI CLC AAUoLl dL LILMI KJ\ \s YJ l_ IVHJUdb Q Case Studies • E. Coli Ammonium Transport • Cell Cycle Regulation Synthetic Biology: Trichlorpropane Degradation FUNDP Namur, 9.5.2014 BioDiVinE Toolset SpEcAff File BIOCHEMICAL MODEL KINETIC ERS MATHEMATICAL MODEL ume (unique iden Description [text] 3e [variable/consta 1 A Variable 2 a Variable 3 c Variable ▼ Reaction name Chemical equation Rate law Description 1 r1 A->B+C mass action decomplex ation 2 r2 B->C mass action conversion 3 r3 c->a mass action conversion o input: • textual: internal .bio format - ODEs + LTL property • gui: list of chemical reactions; SBML standard • tasks: • rectangular abstraction • parallel LTL model checking • output: • model checking counterexample • 2D reachability visualization File View ř =ifei multi-affine ODE model • kinetic parameters set w.r.t. literature • internal and external pH conditions considered constant • initial conditions set to intervals: AmtB, AmtB : NH3, AmtB : NHA NH3in NHAin NH3ex, NHAex (0,1 • 10"5} (1 • 10~b, 1.1 • 10"b) (2 • 10~ö,2.1 • 10~b> (0,1 • 10"5} 9 abstraction - number of discrete concentration levels considered: AmtB AmtB : NH3 AmtB : NHA NH3in NHAin 7 9 3 8 26 FUNDP Namur, 9.5.2014 109/139 E. Coli Ammonium Transport: Model Settings Settings • mass action kinetics ^> multi-affine ODE model • abstraction - number of discrete concentration levels considered: AmtB AmtB : NH3 AmtB : NHA NH3in NHAin 7_9_3_8 26 • initial conditions set to impose low external ammonium conditions Experiments • find the maximal set of parameter values for the given uknown parameter ensuring the maximal reachable level of internal A//-/3 is 1.1 • 106 mol 9 the employed LTL property: G(NH3in < 1.1 • 106) FUNDP Namur, 9.5.2014 E. Coli Ammonium Transport: Experiments params. intervals of validity time k* (1 • 10"12,2.7 • 106) 30 s (5.2 • 106,1 - 1012) 22 s h (1 • 10"12,3.3 • 106) 33 s kg (1 • 10"12,2.7 • 106) 20 s ^1,6,10 see the paper 19 min J. Barnat, L. Brim, A. Krejci, D. Safranek, A. Streck, M. Vejnar, and T. Vejpustek. "On Parameter Synthesis by Parallel Model Checking". IEEE/ACM Transactions on Computational Biology and Bioinformatics. May-June 2012;9(3):693-705 FUNDP Namur, 9.5.2014 111/139 Genetic Regulation of G\/S Transition • central module controlling G\/S transition of mammalian cells FUNDP Namur, 9.5.2014 112/139 Genetic Regulation of G\/S Transition M. Swat, A. Kel, and H. Herzel, "Bifurcation analysis of the regulatory modules of the mammalian Gl/S transition," Bioinformatics, vol. 20, no. 10, pp. 1506-1511, 2004. FUNDP Namur, 9.5.2014 113/139 Genetic Regulation of G\/S Transition ^££1 =k1e1(pRB>E2Fl) -ipaD[vlt.B] «\E*tpi\ =kp + k2&2(pRtl, E2F1) - ~/EzFl [E2F1] o bistability w.r.t. setting of 7p/?e parameter in the range [0.01,1] a liveness properties FG[E2F1] > 8 and FG[E2F1] < 3 are employed a many false-positive runs arise due to time-convergent behaviour introduced by abstraction • by determining transient rectangles we were able to find acceptable results FUNDP Namur, 9.5.2014 Case study: Biodegradation of Trichloropropane in E. coli TCP DhaA ■» DCP HheC ■> ECH EchA + CPD HheC ■» GDL EchA ■> GLY DhaA31 HheC EchA d[TCP] _ kx-DhaA-lTCP] dt ~ KmA + [TCP] d[DCP] _ kv DhaA--[TCP] dt ~ KmA + [TCP] d[ECH] _ k2-HheC-[DCP] _ dt ~ Kma + [DCP] d[CPD] _ k3-EchA-[ECH] _ dt ~ Km,3 + [ECH] d[GDL] _ kA-HheC-[CPD] _ dt ~ KmA+[CPD] d[GLY] _ kb-HheC-[GDL] dt ~ Km^ + [GDL] k2- HheC -[DCP] Kma + [DCP] k3-EchA-[ECH] Km,3 + [ECH] kA-HheC-[CPD] KmA + [CPD] kb-HheC-[GDL] Km^ + [GDL] biodegradation of toxic substrate and intermediates synthetic pathway utilising enzymes from two other bacteria Rhodococcus rhodochrous NCI MB 13064; Agrobacterium radiobacter AD1 find optimal enzymes concentration balancing metabolic burden and toxicity FUNDP Namur, 9.5.2014 115/139 Case study: Biodegradation of Trichloropropane in E. coli Desired behaviour: "TCP is finally completely degraded and the concentration of intermediates does not exceed given bounds" Formally: x)U(FG [TCP] < y)), y>2 = (([GLY] < y)U(FG [GLY] > x)), cps = (G [DCP] 1 Target values assigned to regulatory contexts for all nodes make a PARAMETER SET (parameterization). R. Thomas and R. d'Ari, CRC Press 1990. Biological feedback. FUNDP Namur, 9.5.2014 123/139 Dynamics as a State Transition Graph FUNDP Namur, 9.5.2014 124/139 Dynamics as a State Transition Graph R. Thomas and R. d'Ari, CRC Press 1990. Biological feedback. FUNDP Namur, 9.5.2014 125/139 Dynamics as a State Transition Graph FUNDP Namur, 9.5.2014 126/139 Dynamics as a State Transition Graph R. Thomas and R. d'Ari, CRC Press 1990. Biological feedback. FUNDP Namur, 9.5.2014 127/139 m u Number of possible parameterizations of a single node is exponential w.r.t. the node's in-degree. (more precisely w.r.t. the number of regulatory contexts) FUNDP Namur, 9.5.2014 128/139 Model Checking-based Methodology Time series measurements Initial parametrizations Static constraints (Sec. 2.2) Model checking Dynamic constraints (Sec. 2.3) Hypotheses formalized in LTL 1 Ranking parametrizations Length and robustness of time series walks (Sec. 3) Visualization Behaviour maps and expression profiles (Sec. 4) Acceptable parametrizations Optimal parametrizations Experimental Design o a prototype tool chain: Parsybone - https://github.com/sybila/Parsybone.git ParameterFilter - https://github.com/sybila/ParameterFilter.git o now unified and available online at http://tremppi.fi.muni.cz o distributed computation of acceptable parameterizations 9 employing witnesses (counterexamples) to rank obtained parameterizations • visualization of the results (export to Cytoscape) FUNDP Namur, 9.5.2014 129/139 Time-series Measurement as a Dynamic Constraint FUNDP Namur, 9.5.2014 130/139 Time-series Measurement as a Dynamic Constraint FUNDP Namur, 9.5.2014 130/139 Time-series Measurement as a Dynamic Constraint Expression of v4 along red path Time-series measurement v2 v4 v5 1 1 1 1 1 1 0 1 1 0 1 1 2 2 1 Encoded in LTL: a(l) = All v, = 1 (t(2) = A/e{i,2,4} vi = 1A A/e{2,5} v> = 0 = A/e{i,2,5} v' = 1 A A/e{3,4} = 2 ^ = a(l)A(a(l)U(a(2)AF(a(3)))) ^2 i 3 j 4 ;> & y s 9 mil Discrete time monotonicity between 1st and 2nd measurement 1 5 3 * S 6 y 5 311)11 Discrete time FUNDP Namur, 9.5.2014 130/139 Time-series Measurement as a Dynamic Constraint Expression of v4 along red path Time-series measurement v1 v2 v3 v4 v5 ra 1 1 1 ■ 1 0 1 1 0 1 1 2 2 1 Encoded in LTL: = A/=2 Vi = 1 = 0 = A/6{1,2,5} W = 1 A A/e{3,4} W = 2 ^ = a(l)A(a(l)U(a(2)AF(a(3)))) 0) &2 3* i s 34 § 6 ? a 5iuri Discrete time monotonicity betweerMst and 2n,d measurement 1 3 3 4 & 5 y 5 5 inn Discrete time FUNDP Namur, 9.5.2014 Model Checking on Coloured Graphs • explicit representation of indexed parameter sets (ordered bit vectors) • parameter space split to exclusive blocks equal to size of integer type • each block contains "close" parameter sets • data-parallel distribution: blocks evenly distributed over the cluster 1. V :— color(s) 3. color(t) :— color(t)\V Pi-i Pi Pi+i 2.V: = Pkcolor(s t) FUNDP Namur, 9.5.2014 131/139 Parameterization Ranking: Length Cost • theoretically infinitely many time-series walks • fix a dynamic constraint and focus on compatible shortest walks • penalize unnecessarily higher energy cost • avoid complex model realizations of the constraint 9 assign each parameterization its length cost - the length of a shortest time-series walk 9 consider parameterizations with minimum length cost FUNDP Namur, 9.5.2014 132/139 Parameterization Ranking: Robustness • non-deterministic dynamics caused by asynchronicity • how can we interpret walks with less options to walk off the "optimal path" and miss the expected final state of the time-series? • the property of the model, but... • another classification of parameterizations • local robustness: property Of a State - "umber of valid successors r r J out degree • global robustness: property of a walk - product of local robustness over all states of the walk • model robustness: property of a parameterization - average of global robustness over all time-series walks FUNDP Namur, 9.5.2014 133/139 Parameterization Ranking: Robustness • non-deterministic dynamics caused by asynchronicity • how can we interpret walks with less options to walk off the "optimal path" and miss the expected final state of the time-series? • the property of the model, but... • another classification of parameterizations • local robustness - approximated: 1 Prob(x) =---— out-degree (x J • global robustness: property of a walk - product of local robustness over all states of the walk • model robustness: property of a parameterization - average of global robustness over all time-series walks FUNDP Namur, 9.5.2014 133/139 Parameterization Ranking: Overall Procedure INPUT: regulatory network, initial parameter space, static and dynamic constraints OUTPUT: subset of the initial parameter space containing optimal parameterizations O Remove parametrizations violating static constraints O Compute parameterizations acceptable by dynamic constraints O Select parametrizations with minimal length cost O Select parametrizations with maximal robustness FUNDP Namur, 9.5.2014 134/139 Visualising Results Behaviour Maps and Expression Profiles FUNDP Namur, 9.5.2014 135/139 Case Studies Bacteriophage A1 mon+ v - cl ~+ cro 1 3/ i 1 ell -« [Thieffry et al. 1995] Rat neural system' C2 [Wahde et al. 2001] Init. Parameter Space 6.9 • 109 2.6 • 105 Static Constraints 8.2 • 104 162 Dynamic Constraints 537 108 Length Cost (min) 28 (length 9) 108 (length 5) Robustness (max) 3 (9.7%) 4 (75%) CMSB 2012 Proceedings 2 Fl ML) Technical Report FUNDP Namur, 9.5.2014 136/139 Rat Neural System: Inferring New Hypothesis [Wan 1998, Wahde 2001] Shortest paths Maximally robust paths (10,0,00 Predicted Hypothesis Genes in cluster 4 express before the cluster 1 expression starts to degrade. FUNDP Namur, 9.5.2014 137/139 9 CTL coloured model checking - Pithya Tool https://doi.org/10.1007/978-3-319-63387-9.29 http://pithya.ics.muni.cz • discrete bifurcation analysis checking https://doi.org/10.1016/j.entcs.2015.06.008 • stochastic modelling and parameter synthesis https://doi.org/10.1007/978-3-642-39799-8.7 o robustness analysis https://doi.org/10.1371/j ournal.pone.0094553 • monitoring https://doi.org/10.1016/j.ic.2014.01.012 FUNDP Namur, 9.5.2014 138/139 Thank You for your attention. FUNDP Namur, 9.5.2014 139/139