Nonlinear real arithmetic and -satisfiability Paolo Zuliani School - PowerPoint PPT Presentation

Nonlinear real arithmetic and δ -satisfiability Paolo Zuliani School of Computing Science Newcastle University, UK ( Slides courtesy of Sicun Gao, MIT) 1 / 26

Introduction ◮ We use hybrid systems for modelling and verifying biological system models ◮ prostate cancer therapy ◮ psoriasis UVB treatment ◮ Hybrid systems combine continuous dynamics with discrete state changes 2 / 26

Why Nonlinear Real Arithmetic and Hybrid Systems? (I) A prostate cancer model 1 dx � α x β x � 1 − z � � dt = 1 + e ( k 1 − z ) k 2 − 1 + e ( z − k 3 ) k 4 − m 1 − c 1 x + c 2 z 0 dy � 1 − z � � � z � � dt = m 1 x + α y 1 − d 0 − β y y z 0 z 0 dz dt = − z γ − c 3 v = x + y ◮ v - prostate specific antigen (PSA) ◮ x - hormone sensitive cells (HSCs) ◮ y - castration resistant cells (CRCs) ◮ z - androgen 1A.M. Ideta, G. Tanaka, T. Takeuchi, K. Aihara: A mathematical model of intermittent androgen suppression for prostate cancer. Journal of Nonlinear Science , 18(6), 593–614 (2008) 3 / 26

𝑤 1 = 𝑤 𝑛𝑏𝑦 𝑇 1 ≥ 𝑇 2 + 𝑤 01 ∗ 𝑢 𝑡𝑏𝑔𝑓 𝑤 1 = 0 𝑢 = 𝑢 𝑠𝑓𝑏𝑑𝑢 𝑒𝑇 1 𝑒𝑇 1 𝑒𝑇 1 𝑒𝑇 1 𝑒𝑇 2 𝑒𝑢 = 𝑤 01 𝑒𝑢 = 𝑤 01 + 𝑏 𝑏1 ∗ 𝑢 𝑒𝑢 = 𝑤 01 + 𝑏 𝑒1 ∗ 𝑢 𝑒𝑢 = 𝑤 01 + 𝑏 𝑒1 ∗ 𝑢 𝑒𝑢 = 𝑤 02 + 𝑏 𝑒2 ∗ 𝑢 𝑒𝑇 2 𝑒𝑇 2 𝑒𝑇 2 𝑒𝑇 2 𝑒𝑤 2 𝑒𝑢 = 𝑤 02 𝑒𝑢 = 𝑤 02 𝑒𝑢 = 𝑏 𝑒2 𝑒𝑢 = 𝑤 02 𝑒𝑢 = 𝑤 02 + 𝑏 𝑒2 ∗ 𝑢 𝑒𝑤 1 𝑒𝑤 1 𝑒𝑤 1 𝑒𝑢 = 𝑏 𝑏1 𝑒𝑢 = 𝑏 𝑒1 𝑒𝑢 = 𝑏 𝑒1 𝑒𝑤 2 𝑒𝑢 = 𝑏 𝑒2 Why Nonlinear Real Arithmetic and Hybrid Systems? (I) Intermittent androgen deprivation therapy on-therapy 𝑒𝑦 𝛽 𝑦 1 + 𝑓 𝑨−𝑙 3 𝑙 4 − 𝑛 1 1 − 𝑨 𝛾 𝑦 𝑒𝑢 = 1 + 𝑓 𝑙 1 −𝑨 𝑙 2 − − 𝑑 1 𝑦 + 𝑑 2 𝑨 0 𝑒𝑧 𝑒𝑢 = 𝑛 1 1 − 𝑨 𝑦 + 𝛽 𝑧 1 − 𝑒 0 𝑨 − 𝛾 𝑧 𝑨 0 𝑨 0 𝑒𝑨 𝑒𝑢 = −𝑨𝛿 + 𝑑 3 𝑦 + 𝑧 ≤ 𝑠 0 𝑦 + 𝑧 ≥ 𝑠 1 off-therapy 𝑒𝑦 𝛽 𝑦 1 + 𝑓 𝑨−𝑙 3 𝑙 4 − 𝑛 1 1 − 𝑨 𝛾 𝑦 𝑒𝑢 = 1 + 𝑓 𝑙 1 −𝑨 𝑙 2 − − 𝑑 1 𝑦 + 𝑑 2 𝑨 0 𝑒𝑧 𝑒𝑢 = 𝑛 1 1 − 𝑨 𝑦 + 𝛽 𝑧 1 − 𝑒 0 𝑨 − 𝛾 𝑧 𝑨 0 𝑨 0 𝑒𝑨 𝑒𝑢 = 𝑨 0 − 𝑨 𝛿 + 𝑑 3 4 / 26

Why Nonlinear Real Arithmetic and Hybrid Systems? (II) A model of psoriasis development and UVB treatment 2 ω (1 − SC + λ SCd ) SC dSC k 1 s ω SCmax = γ 1 − β 1 In A SC − 1 + ( ω − 1)( TA + TAd 1 + ( ω − 1)( TA + TAd dt ) n ) n SC + k 1 TA Pta , h Pta , h dTA k 1 a , s ω SC 2 k 1 s ω = + 1 + ( ω − 1)( TA + TAd 1 + ( ω − 1)( TA + TAd dt ) n + γ 2 GA − β 2 In A TA − k 2 s TA − k 1 TA ) n Pta , h Pta , h dGA = ( k 2 a , s + 2 k 2 s ) TA − k 2 GA − k 3 GA − β 3 GA dt k p SC 2 dSC d SC + SC d d = γ 1 d (1 − + k 1 d TA d ) SC d − β 1 d In A SC d − k 1 sd SC d − a + SC 2 k 2 dt SC max , t d dTA d = k 1 a , sd SC d + 2 k 1 sd SC d + γ 2 d TA d + k 2 d GA d − β 2 d In A TA d − k 2 sd TA d − k 1 d TA d dt dGA d = ( k 2 a , sd + 2 k 2 sd ) TA d − k 2 d GA d − k 3 d GA d − β 3 d GA d dt ◮ Therapy episode: 48 hours of irradiation + 8 hours of rest 2H. Zhang, W. Hou, L. Henrot, S. Schnebert, M. Dumas, C. Heus` ele, and J. Yang. Modelling epidermis homoeostasis and psoriasis pathogenesis. Journal of The Royal Society Interface , 12(103), 2015. 5 / 26

Why Nonlinear Real Arithmetic and Hybrid Systems? (II) A model of psoriasis development and UVB treatment 2 ω (1 − SC + λ SCd ) SC dSC k 1 s ω SCmax = γ 1 − β 1 In A SC − 1 + ( ω − 1)( TA + TAd 1 + ( ω − 1)( TA + TAd dt ) n ) n SC + k 1 TA Pta , h Pta , h dTA k 1 a , s ω SC 2 k 1 s ω = + 1 + ( ω − 1)( TA + TAd 1 + ( ω − 1)( TA + TAd dt ) n ) n + γ 2 GA − β 2 In A TA − k 2 s TA − k 1 TA Pta , h Pta , h dGA = ( k 2 a , s + 2 k 2 s ) TA − k 2 GA − k 3 GA − β 3 GA dt k p SC 2 dSC d SC + SC d d = γ 1 d (1 − SC d − β 1 d In A SC d − k 1 sd SC d − + k 1 d TA d ) k 2 a + SC 2 dt SC max , t d dTA d = k 1 a , sd SC d + 2 k 1 sd SC d + γ 2 d TA d + k 2 d GA d − β 2 d In A TA d − k 2 sd TA d − k 1 d TA d dt dGA d = ( k 2 a , sd + 2 k 2 sd ) TA d − k 2 d GA d − k 3 d GA d − β 3 d GA d dt ◮ Therapy episode: 48 hours of irradiation + 8 hours of rest ◮ Therapy episode = multiply β 1 and β 2 by a constant In A 2H. Zhang, W. Hou, L. Henrot, S. Schnebert, M. Dumas, C. Heus` ele, and J. Yang. Modelling epidermis homoeostasis and psoriasis pathogenesis. Journal of The Royal Society Interface , 12(103), 2015. 5 / 26

Real-World Applications ◮ Psoriasis Stratification to Optimise Relevant Therapy (PSORT) ◮ Large ( ∼ £ 5m) ◮ Primarily biomarkers discovery ◮ We use computational modelling for understanding psoriasis’ mechanisms 6 / 26

Real-World Applications ◮ Psoriasis Stratification to Optimise Relevant Therapy (PSORT) ◮ Large ( ∼ £ 5m) ◮ Primarily biomarkers discovery ◮ We use computational modelling for understanding psoriasis’ mechanisms ◮ Personalised ultraviolet B treatment of psoriasis through biomarker integration with computational modelling of psoriatic plaque resolution ◮ Starts February 2017 ◮ PIs: P.Z. and Nick Reynolds (Institute of Cellular Medicine) ◮ Computational modelling to inform UVB therapies used in the clinic — real impact on people’s health! 6 / 26

Bounded Reachability ◮ Reachability is a key property in verification, also for hybrid systems ◮ Reachability is undecidable even for linear hybrid systems (Alur, Courcoubetis, Henzinger, Ho. 1993) ◮ [ Bounded Reachability] Does the hybrid system reach a goal state within a finite time and number of (discrete) steps? 7 / 26

Bounded Reachability ◮ Reachability is a key property in verification, also for hybrid systems ◮ Reachability is undecidable even for linear hybrid systems (Alur, Courcoubetis, Henzinger, Ho. 1993) ◮ [ Bounded Reachability] Does the hybrid system reach a goal state within a finite time and number of (discrete) steps? ◮ “ Can a 5-episode UVB therapy remit psoriasis for a year ?” 7 / 26

Bounded Reachability ◮ Reachability is a key property in verification, also for hybrid systems ◮ Reachability is undecidable even for linear hybrid systems (Alur, Courcoubetis, Henzinger, Ho. 1993) ◮ [ Bounded Reachability] Does the hybrid system reach a goal state within a finite time and number of (discrete) steps? ◮ “ Can a 5-episode UVB therapy remit psoriasis for a year ?” ◮ Reasoning about nonlinear real arithmetic is hard . . . 7 / 26

Type 2 Computability Turning machines operate on finite strings, i.e. , integers, which cannot capture real-valued functions. ◮ Real numbers can be encoded on infinite tapes. ◮ Real numbers are functions over integers. ◮ Real functions can be computed by machines that take infinite tapes as inputs, and output infinite tapes encoding the values. Definition (Name of a real number) A real number a can be encoded by an infinite sequence of rationals γ a : N → Q such that ∀ i ∈ N | a − γ a ( i ) | < 2 − i . 8 / 26

Type 2 Computability A function f ( x ) = y is computable if any name of x can be algorithmically mapped to a name of y y 1 } ... . . k input tapes . y k ... } ... ... . . work tapes . M ... ... y output tape ... f M ( y 1 , . . . , y k ) = y Writing on any finite segment of the output tape takes finite time. 9 / 26

Type 2 Computability ◮ Type 2 computability implies continuity ◮ “Numerically computable” roughly means Type 2 computable ◮ Approximation up to arbitrary numerical precisions Ker-I Ko. Complexity Theory of Real Functions . 1991. 10 / 26

Facts Type 2 Computable: ◮ polynomials, sin, exp, . . . ◮ numerically feasible ODEs, PDEs, . . . Type 2 Complexity: ◮ sin, exp, etc. are in P [0 , 1] ◮ Lipschitz-continuous ODEs are in PSPACE [0 , 1] ; in fact, can be PSPACE [0 , 1] -complete (Kawamura, CCC 2009). See Ko’s book for many more results . . . 11 / 26

L R F -Formulas (Gao, Avigad, and Clarke. LICS 2012) Let F be the class of all Type 2 computable real functions. Definition ( L R F -Formulas) First-order language over � >, F� : t := x | f ( t ( � x )) ϕ := t ( � x ) > 0 | ¬ ϕ | ϕ ∨ ϕ | ∃ x i ϕ | ∀ x i ϕ Example Let dx / dt = f ( x ) be an n-dimensional dynamical system. Lyapunov stability is expressed as: � t � � ∀ ε ∃ δ ∀ t ∀ x 0 ∀ x t . || x 0 || < δ ∧ x t = x 0 + f ( s ) ds → || x t || < ε 0 12 / 26

Hybrid Automata A hybrid automaton is a tuple y ) : q , q ′ ∈ Q } , H = � X , Q , { flow q ( � x , � y , t ) : q ∈ Q } , { jump q → q ′ ( � x , � { inv q ( � x ) : q ∈ Q } , { init q ( � x ) : q ∈ Q }� ◮ X ⊆ R n for some n ∈ N ◮ Q = { q 1 , ..., q m } is a finite set of modes ◮ Other components are finite sets of quantifier-free L R F -formulas. 13 / 26

Example: Nonlinear Bouncing Ball ◮ X = R 2 and Q = { q u , q d } . ◮ flow q d ( x 0 , v 0 , x t , v t , t ), dynamics in the falling phase: � t � t g (1 + β v ( s ) 2 ) ds ) ( x t = x 0 + v ( s ) ds ) ∧ ( v t = v 0 + 0 0 ◮ jump q u → q d ( x , v , x ′ , v ′ ): ( v = 0 ∧ x ′ = x ∧ v ′ = v ) ◮ inv q d : ( x > = 0 ∧ v > = 0). ◮ init q d : ( x = 10 ∧ v = 0). 14 / 26