Grégory Batt Ph.D. defense -- under supervision of Hidde de Jong, Helix research group

1. Validation of Qualitative Models of Genetic Regulatory NetworksA Method Based on Formal Verification Techniques Grégory Batt Ph.D. defense -- under supervision of Hidde de Jong, Helix research group INRIA Rhône-Alpes -- Ecole doctorale Mathématiques, Sciences et technologies de l’information, Informatique Université Joseph Fourier

2. Heat shock Nutritional stress Cold shock Osmotic stress … Stress response in Escherichia coli • Bacteria capable of adapting to a variety of changing environmental conditions • Stress response in E. coli has been much studied Model for understanding adaptation of pathogenic bacteria to their host

3. Nutritional stress response in E. coli • Response of E. coli to nutritional stress conditions: transition from exponential phase to stationary phase Important developmental decision: profound changes of morphology, metabolism, gene expression,... log (pop. size) > 4 h time

4. fis P gyrAB P cya P1-P’1 P2 FIS GyrAB CYA DNA supercoiling cAMP•CRP Signal (carbon starvation) TopA CRP stable RNAs topA P1-P4 crp P1 P2 protein rrn P1 P2 gene promoter Network controlling stress response • Response of E. coli to nutritional stress conditions controlled by genetic regulatory network Despite abundant knowledge on network components, no global view of functioning of network available

5. Modeling and simulation • Genetic regulatory network controlling E. coli stress response is large and complex • Modeling and simulation indispensable for dynamical analysis of genetic regulatory networks Systematic prediction of possible network behaviors • Current constraints on modeling and simulation: • knowledge on molecular mechanisms rare • quantitative information on kinetic parameters and molecular concentrations absent • Qualitative methods developed for analysis of genetic networks using coarse-grained models

6. . x = f (x) consistency? experimental data model network predictions Model validation • Available information on structure of network controlling E. coli stress response is incomplete Model is working hypothesis and needs to be tested • Model validation is prerequisite for use of model as predictive and explanatory tool Check consistency between model predictions and experimental data

7. Model validation • Available information on structure of network controlling E. coli stress response is incomplete Model is working hypothesis and needs to be tested • Model validation is prerequisite for use of model as predictive and explanatory tool Check consistency between model predictions and experimental data • Current constraints on model validation: • available experimental data essentially qualitative in nature • model validation must be automatic and efficient

8. Objectives and approach of thesis • Objective of thesis: Development of automated and efficient method for testing whether predictions from qualitative models of genetic regulatory networks are consistent with experimental data on dynamics of system • Approach based on formal verification of hybrid systems • qualitative analysis of piecewise-linear models of genetic networks • model checking for testing consistency between predictions and data • Expected contributions: • scalable method with sound theoretical basis • implementation of method in user-friendly computer tool • applications to validation of models of networks of biological interest

9. Overview • Introduction • Method for model validation • Piecewise-linear (PL) differential equation models • Symbolic analysis using qualitative abstraction • Verification of properties by means model-checking techniques • Genetic Network Analyzer 6.0 • Validation of model of nutritional stress response in E. coli • Discussion and conclusions

10. Overview • Introduction • Method for model validation • Piecewise-linear (PL) differential equation models • Symbolic analysis using qualitative abstraction • Verification of properties by means model-checking techniques • Genetic Network Analyzer 6.0 • Validation of model of nutritional stress response in E. coli • Discussion and conclusions

11. s-(x, θ) 1 0 x  A x : protein concentration  : threshold concentration . . B  ,  : rate constants xbbs-(xa, a1) – bxb xaas-(xa, a2) s-(xb, b ) – axa b a PL differential equation models • Genetic networks modeled by class of differential equations using step functions to describe switch-like regulatory interactions • Hybrid, piecewise-linear (PL) models of genetic regulatory networks Glass and Kauffman, J. Theor. Biol., 73

12. maxb M15 M14 M13 M12 M11 M8 M6 M7 M10 M9 b M1 M2 M3 M4 M5 maxb maxb 0 maxa a1 a2 b b 0 0 maxa maxa a1 a1 a2 a2 Qualitative analysis of network dynamics . x = h (x), x  \ • Analysis of the dynamics in phase space: • Partition of phase space into mode domains maxb b 0 maxa a1 a2

13. kb/gb maxb b ka/ga . . . . xaa – axa xbb – bxb xbbs-(xa, a1) – bxb xaas-(xa, a2) s-(xb, b ) – axa 0 < qa1 < qa2 < a/a < maxa 0 maxa a1 a2 0 a – axa 0 < qb < b/b< maxb 0 b – bxb Qualitative analysis of network dynamics . x = h (x), x  \ • Analysis of the dynamics in phase space: maxb b M1 0 maxa a1 a2

14. kb/gb maxb b 0 < qa1 < qa2 < a/a < maxa 0 maxa a1 a2 0 < qb < b/b< maxb Qualitative analysis of network dynamics . x = h (x), x  \ • Analysis of the dynamics in phase space: maxb . M11 xa– axa . xbb – bxb b 0 maxa a1 a2

15. kb/gb maxb M1 M2 M3 b ka/ga 0 maxa a1 a2 Qualitative analysis of network dynamics . x = h (x), x  \ • Analysis of the dynamics in phase space: maxb b 0 maxa a1 a2

16. maxb maxb b b 0 0 maxa maxa a1 a1 a2 a2 Qualitative analysis of network dynamics . x = h (x), x  \ • Analysis of the dynamics in phase space: • Extension of PL differential equations to differential inclusions using Filippov approach: maxb maxb kb/gb b b M3 M4 M1 M2 M3 M5 0 0 maxa a1 maxa ka/ga a1 a2 ka/ga a2 . x H (x), x  • Gouzé and Sari, Dyn. Syst., 02

17. maxb M15 M14 M13 M12 M11 M8 M6 M7 M10 M9 b M1 M2 M3 M4 M5 maxb maxb 0 maxa a1 a2 b b 0 0 maxa maxa a1 a1 a2 a2 Qualitative analysis of network dynamics . • Analysis of the dynamics in phase space: • In every mode domain M, the system either converges monotonically towards focal set, or instantaneously traverses M x H (x), x  maxb b 0 maxa a1 a2 de Jong et al., Bull. Math. Biol., 04 • Gouzé and Sari, Dyn. Syst., 02

18. maxb maxb M15 M14 M13 M12 M11 M8 M6 M7 M10 M9 b b M1 M2 M3 M4 M5 maxb maxb 0 0 maxa maxa a1 a1 a2 a2 b b . . x M11: xa < 0, xb? 0 0 maxa maxa a1 a1 a2 a2 Problem for model validation • Partition does not preserve unicity of derivative sign Predictions not adapted to comparison with available experimental data: temporal evolution of direction of change of protein concentrations

19. maxb b maxb maxb D21 D24 D20 kb/gb D18 D19 D23 D27 D26 D25 0 D17 D16 D22 maxa a1 a2 b b D14 D13 D15 D11 D10 D12 . . x D17: xa < 0, xb > 0 D1 D3 D5 D7 D9 D2 D4 D6 D8 0 0 maxa maxa a1 a1 a2 a2 Qualitative analysis of network dynamics • Finer partition of phase space: flow domains • In every domain D, the system either converges monotonically towards focal set, or instantaneously traverses D • In every domain D, derivative signs are identical everywhere maxb b 0 maxa a1 a2

20. Continuous transition system • PL system, = (,,H), associated with continuous PL transition system,-TS = (,→,╞), where • continuous phase space

21. : transition from x to x’ iff a solution reaches x’ from x and xand x’ in same or in adjacent domain maxb maxb x5 kb/gb x4 x3 b b x1 x2 0 0 a1 a1 a2 a2 maxa maxa Continuous transition system • PL system,  = (,,H), associated with continuous PL transition system,-TS = (,→,╞), where • continuous phase space • →transition relation x1 → x2, x1→ x3, x3→ x4 x2→ x3,

22. . . . . x4╞xb> 0, x4╞xa< 0, x1╞xb> 0, x1╞xa> 0, : describes derivative sign of solutions at x Continuous transition system • PL system,  = (,,H), associated with continuous PL transition system,-TS = (,→,╞), where • continuous phase space • →transition relation • ╞satisfaction relation • and -TShave equivalent reachability properties maxb maxb x5 kb/gb x4 x3 b b x1 x2 0 0 a1 a1 a2 a2 maxa maxa

23. maxb maxb D21 D24 D20 D18 D19 D23 D27 D26 D25 kb/gb D17 D16 D22 D14 D13 D15 D11 D10 D12 b b D1 D3 D5 D7 D9 D2 D4 D6 D8 D1 D ; 0 0 a1 a1 a2 a2 maxa maxa Discrete abstraction • Qualitative PL transition system,-QTS = (D, →,╞), where • D finite set of domains :D= {D1, …, D27}

24. : transition from D to D’ iff there exist xD, x’D’ such that x→ x’ maxb maxb x5 kb/gb x4 D17 x3 D17 b b D11 x1 x1 x2 D11 D1 D1 D ; D1 →~D1, D1 →~D11, D11 →~D17, D1 D1 0 0 a1 a1 a2 a2 maxa maxa Discrete abstraction • Qualitative PL transition system,-QTS = (D, →,╞), where • D finite set of domains • → quotient transition relation

25. maxb maxb x5 kb/gb x4 x3 D17 b b D11 x1 x1 x2 D1 0 0 a1 a1 a2 a2 maxa maxa . . . D1 D ; D1 →~D1, D1 →~D11, D11 →~D17, D1╞ xa>0, D1╞xb>0,D4╞ xa < 0 Discrete abstraction • Qualitative PL transition system,-QTS = (D,→,╞), where • D finite set of domains • → quotient transition relation • ╞ quotient satisfaction relation : D╞piff there exists xD such thatx╞p

26. D21 D24 D20 D18 D18 D19 D23 D27 D26 D25 D21 D24 D20 D17 D17 D16 D22 D25 D27 D26 D19 D23 D14 D13 D15 D11 D10 D12 D18 D11 D17 D22 D16 D1 D3 D5 D7 D9 D2 D4 D6 D8 D13 D11 D10 D15 D12 D14 D1 D1 D9 D3 D5 D7 D2 D4 D6 D8 Discrete abstraction • Qualitative PL transition system,-QTS = (D,→,╞), where • D finite set of domains • → quotient transition relation • ╞quotient satisfaction relation • Quotient transition system -QTS is a simulation of-TS (but not a bisimulation) maxb maxb kb/gb b b 0 0 a1 a1 a2 a2 maxa maxa Alur et al., Proc. IEEE, 00

27. Discrete abstraction • Important properties of -QTS : • -QTS provides finite and qualitative description of the dynamics of system  in phase space • -QTS is aconservative approximation of : every solution of corresponds to a path in -QTS • -QTS is invariant for all parameters , ,and  satisfying a set of inequality constraints • -QTS can be computed symbolically using parameter inequality constraints: qualitative simulation • Use of-QTS to verify dynamical properties of original system  Need for automatic and efficient method to verify properties of -QTS Batt et al., HSCC, 05

28. Model-checking approach • Model checking is automated technique for verifying that discrete transition system satisfies certain temporal properties • Computation tree logic model-checking framework: • set of atomic propositionsAP • discrete transition system is Kripke structure KS =  S, R, L , where Sset of states, Rtransition relation, Llabeling function over AP • temporal properties expressed in Computation Tree Logic (CTL) p, ¬f1, f1f2, f1f2, f1→f2, EXf1, AXf1, EFf1, AFf1, EGf1, AGf1, Ef1Uf2, Af1Uf2, where pAP and f1, f2 CTL formulas • Computer tools are available to perform efficient and reliable model checking (e.g., NuSMV, SPIN, CADP)

29. . . . xa . . . . xa > 0 xb > 0 xb > 0 xa < 0 0 time xb . . There Exists a Future state wherexa > 0and xb > 0 and from that state, thereExists a Future state wherexa < 0and xb > 0 0 . . time . . . . EF(xa > 0 xb > 0EF(xa < 0 xb > 0) ) Validation using model checking • Atomic propositions AP ={xa = 0, xa <qa1, ... , xb < maxb, xa< 0, xa= 0, ... , xb> 0} • Observed property expressed in CTL

30. D21 D24 D20 D25 D27 D26 D19 D23 xa D18 . . . . D17 D22 xa < 0 xb > 0 xb > 0 xa > 0 D16 0 time xb D13 D11 D10 D15 D12 D14 0 time D1 D9 D5 D7 D3 D6 D2 D4 . . . . D8 EF(xa > 0 xb > 0EF(xa < 0 xb > 0) ) Validation using model checking • Discrete transition system computed using qualitative simulation • Use of model checkers to check consistency between experimental data and predictions • Fairness constraints used to exclude spurious behaviors Consistency? Yes Batt et al., IJCAI, 05

31. Overview • Introduction • Method for model validation • Piecewise-linear (PL) differential equation models • Symbolic analysis using qualitative abstraction • Verification of properties by means model-checking techniques • Genetic Network Analyzer 6.0 • Validation of model of nutritional stress response in E. coli • Discussion and conclusions

32. Integration into environment for explorative genomics at Genostar SA Genetic Network Analyzer • Model validation approach implemented in version 6.0 of GNA, freely available for academic research Batt et al., Bioinformatics, 05

33. structure into packages class diagram of kernel Genetic Network Analyzer • GNA implemented in Java 1.4 • > 17000 lines of code in 6 packages • 35% of lines modified with respect to version 5.5 (up to 60% in kernel)

34. Genetic Network Analyzer • Rules for symbolic computation of refined partition and corresponding transition relation and domain properties Tailored algorithms and implementation favor upscalability • Export functionalities to model checkers (NuSMV, CADP)

35. Overview • Introduction • Method for model validation • Piecewise-linear (PL) differential equation models • Symbolic analysis using qualitative abstraction • Verification of properties by means model-checking techniques • Genetic Network Analyzer 6.0 • Validation of model of nutritional stress response in E. coli • Discussion and conclusions

36. Nutritional stress response in E. coli • Entry into stationary phase is an important developmental decision exponential phase stationary phase ? signal of nutritional deprivation How does lack of nutrients induce decision to stop growth?

37. P gyrAB P fis P1-P’1 P2 cya Fis GyrAB CYA Signal (carbon starvation) Supercoiling cAMP•CRP TopA CRP stable RNAs P1-P4 topA crp P1 P2 rrn P1 P2 Model of nutritional stress response • Carbon starvation network modeled by PL model 7 PL differential equations, 40 parameters and 54 inequality constraints Ropers et al.,Biosystems, in press How does response emerge from network of interactions?

38. “Fis concentration decreases and becomes steady in stationary phase” Ali Azam et al., J. Bacteriol., 99 . . EF(xfis < 0EF(xfis = 0 xrrn < qrrn) ) Validation of stress response model • Qualitative simulation of carbon starvation: • 66 reachable domains (< 1s.) • single attractor domain (asymptotically stable equilibrium point) • Experimental data on Fis: CTL formulation: Model checking with NuSMV: property true (< 1s.)

39. . AG(xcrp > q3crp  xcya > q3cya xs > qs → EF xcya < 0) . . EF( (xgyrAB < 0 xtopA > 0) xrrn < qrrn) Validation of stress response model • Other properties: • “cya transcription is negatively regulated by the complex cAMP-CRP” • “DNA supercoiling decreases during transition to stationary phase” • Inconsistency between observation and prediction calls for model revision or model extension Nutritional stress response model extended with global regulator RpoS Kawamukai et al., J. Bacteriol., 85 True (<1s) Balke and Gralla, J. Bacteriol., 87 False (<1s)

40. Novel prediction of stress response model • Qualitative simulation of carbon upshift response: • 1143 reachable domains (< 2s) • several strongly connected components • Are some strongly connected components attractors? • Attractor corresponds to damped oscillations towards stable equilibrium point: unexpected prediction • Experimental verification of model predictions Time-series measurements of protein concentrations in parallel and at high sampling rate using gene reporter system AG(statesInSCCi→ AG statesInSCCi) True (<1s, i=3) Grognard et al., in preparation

41. Overview • Introduction • Method for model validation • Piecewise-linear (PL) differential equation models • Symbolic analysis using qualitative abstraction • Verification of properties by means model-checking techniques • Genetic Network Analyzer 6.0 • Validation of model of nutritional stress response in E. coli • Discussion and conclusions

42. Summary • Development of automated and efficient method for testing whether predictions from qualitative models of genetic regulatory networks are consistent with experimental data on system dynamics • Use of discrete abstraction that yields predictions well-adapted to comparison with available experimental data • Combination of tailored symbolic analysis and model checking for verification of dynamical properties of hybrid models of large and complex networks • Biological relevance demonstrated on validation of models of networks of biological interest Batt et al., HSCC, 05 Batt et al., IJCAI, 05 Batt et al., Bioinformatics, 05

43. Discussion • Discrete abstractions used for analysis of continuous and hybrid models • symbolic reachability analysis of hybrid automata models • more precise analysis of system dynamics • need for complex decision procedures • no treatment of discontinuities in vector field • qualitative simulation using qualitative differential equations • more general class of model • methods are not scalable • Model checking used for analysis of discrete models • verification of properties of logical models • intuitive connection between underlying continuous dynamics and discrete representation • no explicit representation of dynamical phenomena at threshold concentrations Ghosh and Tomlin, Systems Biology, 04 Heidtke and Schulze-Kremer, Bioinformatics, 98 Bernot et al., J. Theor. Biol., 04

44. Perspectives • Further integration of model-checking task into GNA Property specification, verification, interpretation of diagnostics • Exploitation of advanced model-checking techniques Partial order reduction, graph minimization, modular model checking, ... • Extensions of model validation • model inference: complete partially-specified models • model revision: modify inconsistent models • network design: find model satisfying set of design constraints

45. Thanks for your attention!