Dissecting Multiscale Systems Peter Szmolyan Technische Universit - - PowerPoint PPT Presentation

Dissecting Multiscale Systems

Peter Szmolyan

Technische Universit¨ at Wien, Austria

Introduction

• 1. Overview: time scales, slow-fast

dynamics, model reduction

• 2. Olsen model: complicated dynamics
• 3. Cell cycle model: Mitotic Oscillator
• 4. Extensions to other regulatory or

signaling networks

• 5. Conclusion and Outlook

• 1. Biological processes occur on widely different

time scales and show slow-fast dynamics Rhythm Period Neural rhythms 0.01 - 1 s Cardiac rhythm 1 s Calzium-oszillations 1 s – min Biochemical oscillations 1 min – 20 min Cell cycle , signaling processes 10 min – 24 h Hormonal rhythms 10 min – 24 h Circadian rhythm 24 h Ovarian cycle 28 days periodic phenomena,

• A. Goldbeter (1996)

Chemical, metabolic, and gene regulatory networks are often modeled as systems of ODEs many variables many parameters very different orders of magnitude numerical simulation not too difficult mathematical analysis challenging simple dynamics complicated dynamics many different time scales, slow-fast dynamics

Concept of slow manifolds is widely used

dynamical system: ˙ X = F(X), X 2 M all solutions rapidly approach lower-dimensional invariant slow manifold S ⇢ M, “fast modes are slaved to slow modes”, mathematical justification: center-, slow-, inertial-manifolds

Often several - hard to identify - slow manifolds exist Difficulties: several branches, different orders of magnitudes different dimensions, different types of stability, saddle type slow manifolds, bifurcations as parameters vary intersections? connections by fast dynamics?

Existing methods for dimension reduction need to be combined quasi steady state assumptions ad hoc reductions, lumping,... classical singular perturbation methods: expansions, Tikhonov theory geometric singular perturbation theory: invariant manifolds, Fenichel theory, blow-up method truncations based on tropical methods symbolic methods computational methods

Good scalings are crucial for model simplification and dimension reduction large terms dominate small terms finding a good scaling is nontrivial! what is a good scaling? nonlinear problem ) good scaling depends on position in phase space

• ften there exist several good scalings

Good scalings are crucial for model simplification and dimension reduction large terms dominate small terms finding a good scaling is nontrivial! what is a good scaling? nonlinear problem ) good scaling depends on position in phase space

• ften there exist several good scalings

x0 = x + εx + εx2, x 2 R, ε ⌧ 1 x = O(1) = ) x0 = x + O(ε) x = O(ε1), x = X ε = ) X0 = X + X2 + O(ε)

Most scalings have the form of power-laws x0 = f(x, µ), x 2 Rn, µ 2 Rm variables x, parameters µ reference value x = ξ, µ = m, often ξ = 0, m = 0 xi = ξi + εαi ˜ xi, i = 1, . . . , n µj = mj + εβj ˜ µj, j = 1, . . . , m key issue: find suitable ε ⌧ 1, α 2 Qn, β 2 Qm tools: experience, intuition, Newton polygon, tropical methods, symbolic computations, numerical simulations,...

Regular perturbation problems are easy, singular perturbation problems are interesting A perturbation problem ˙ u = F(u, ε), 0 < ε ⌧ 1 is regular, if it can be analyzed in a single scaling singular, if two or more different scalings are needed.

Singularly perturbed (slow-fast) ODEs in standard form require two scalings ˙ x = f(x, y, ε) ε ˙ y = g(x, y, ε) (1) x slow, y fast, ε ⌧ 1, slow time scale t, transform to fast time scale τ := t/ε x0 = εf(x, y, ε) y0 = g(x, y, ε) (2)

• Syst. (1) and Syst. (2) equivalent for ε > 0

There are two distinct limiting systems for ε = 0 reduced problem ˙ x = f(x, y, 0) 0 = g(x, y, 0) chemistry: quasi steady state approximation layer problem x0 = y0 = g(x, y, 0) chemistry: fast approach to equilibrium critical manifold S := {g(x, y, 0) = 0} reduced problem is a dynamical system on S S is a “manifold” of equilibria for layer problem

Geometric Singular Perturbation Theory (GSPT) based on invariant manifold theory allows to go from ε = 0 to 0 < ε ⌧ 1 Theorem: critical manifold S normally hyperbolic, i.e. ∂g

∂y|S hyperbolic ) S perturbs smoothly to

slow manifold Sε for ε small, ...

• N. Fenichel (1979)

Almost all slow-fast problems involve non-hyperbolic points folds of S self-intersections of S unbounded branches of S can be treated with blow-up method blow-up, well known in algebraic geometry rescaling + compactification simplest blow-up: polar coordinates

Slow-fast decompositions, GSPT + blow-up method have proven to be very useful enzyme kinetics mathematical neuroscience metabolic processes: glycolysis, oxidization of NADH cell cycle signaling pathways

• 2. Olsen model describes oxidization of

Nicotinamide Adenine Dinucleotide (NADH) 2NADH + O2 + 2H+ ! 2NAD+ + 2H2O ˙ A = k7 k9A k3ABY ˙ B = k8 k1BX k3ABY ˙ X = k1BX 2k2X2 + 3k3ABY k4X + k6 ˙ Y = 2k2X2 k3ABY k5Y A oxygen, B NADH, X, Y intermediate products Reaction rates: k1 = 0.16, 0.35, 0.41 k2 = 250, k3 = 0.035, k4 = 20, k5 = 5.35, k6 = 105, k7 = 0.8, k8 = 0.825, k9 = 0.1

The Olsen model has complicated dynamics a) k1 = 0.16, b) k1 = 0.35, c) k1 = 0.41 slow-fast dynamics: a) and b) mixed-mode oscillations or chaotic c) relaxation oscillations

• C. Kuehn and P. Sz., J. Nonlinear Science (2015)

Visualization in phase space shows more details slow dynamics in A, B close to X, Y ⇡ 0, fast dynamics in A, B, X, Y away from X, Y ⇡ 0,

Scaling A, B, X, and Y gives a slow-fast system (o) ˙ a = θ αa aby ˙ b = ν(1 bx2 aby) ε2 ˙ x = bx x2 + 3aby βx + δ ε2 ˙ y = x2 y aby ν ⇡ 101, θ, α, β ⇡ 1, ε2 ⇡ 102, δ ⇡ 105 ε ⌧ 1 = ) a, b slow x, y fast α ⇠ k1, δ ⇠ k6, bifurcation parameters

Finding the scaling is not easy A = k1k5 k3 p2k2k8 a, B = p2k2k8 k1 b X = k8 2k2 x, Y = k8 k5 y T = k1k5 k3k8 p2k2k8 t

• A. Milik (1996)

The Olsen model has a complicated critical manifold ˙ a = θ αa aby ˙ b = ν(1 bx aby) ε2 ˙ x = bx x2 + 3aby βx + δ ε2 ˙ y = x2 y aby For δ = 0 critical manifold S is defined by 0 = bx x2 + 3aby βx 0 = x2 y aby solve for x, y Remark: δ 6= 0 breaks invariance of x = 0, y = 0

S has several branches and singularities (x, y) = (0, 0) attracting b < β, saddle b > β and x = b β 1 2ab(1 + ab), y = (b β)2 (1 2ab)2(1 + ab) attracting b > β, saddle b > β poles at 2ab = 1, transcritical singularities at b = β (non-hyperbolic)

GSPT explains dynamics for x, y = O(1), ε ⌧ 1 reduced flow in x = 0, y = 0 ˙ a = θ αa ˙ b = ν return mechanism? needs x, y large!

Rescaling for x, y large explains return mechanism scaling

x = X ε , y = Y ε2, t = ε2τ

transforms System (o) to (O) a0 = ε2(µ αa) abY b0 = ε2ν εbY νabY εX0 = X2 + ε(b β)X + 3abY + ε2δ Y 0 = X2 Y abY. a, b, Y slow, X fast

System (O) has a new three-dimensional attracting critical manifold C new critical manifold C := {X2 = 3abY } attracting, X > 0 non-hyperbolic at X = 0, Y = 0 reduced problem on C is integrable a0 = abY b0 = νabY Y 0 = (2ab 1)Y

Large amplitude returns occur on critical manifold C defined by 3abY = X2 GSPT ) rigorous results for (X, Y ) 6= (0, 0), ε ⌧ 1 Close to (X, Y ) = (0, 0) matching of System (O) to System (o) by blow-up method.

Blow-up method allows to match System (O) near (X, Y ) = (0, 0) and System (o) with x, y large Blow-up Transformation: include ε as a variable! X = ¯ r¯ x, Y = ¯ r2¯ y, ε = ¯ r¯ ε, (¯ x, ¯ y, ¯ ε) 2 S2, ¯ r 2 R weighted spherical blow-up of (0, 0, 0) in (X, Y, ε) space diffeomorphism for ¯ r > 0 degenerate point (0, 0, 0) is blown-up to S2 ⇥ {0} vector field can be pulled back to ¯ r = 0 and can be desingularized blow-up = rescaling + compacification

Two charts cover the relevant parts of the sphere S2 blow-up: X = ¯ r¯ x, Y = ¯ r2¯ y, ε = ¯ r¯ ε Chart K1 defined by setting ¯ x = 1, covers x > 0 X = r1, Y = r2

1y1, ε = r1ε1

is equivalent to System (O) for r1 > 0 Chart K2 defined by setting ¯ ε = 1, covers ε > 0 X = r2x2, Y = r2

2y2, ε = r2

reproduces System (o): X = εx, Y = ε2y At ¯ r = 0 the system can be desingularized and rigorous matching near (¯ ε, ¯ r) = (0, 0) is possible.

Why does it work? ¯ C gets separated from the bad set (¯ x, ¯ y) = (0, 0); gain of normal hyperbolicity in parts of K1 and K2 we recover Systems (O) and (O2) which can be analysed by GSPT Within K2 other blow-ups have to be used to study transcritical singularities (δ = 0) or folds (δ > 0)

• 3. The mitotic oscillator is a simple model related

to the dynamics of the cell-cycle

Paul Nurse, Lee Hartwell, Tim Hunt Nobel-Prize in medicine (2001)

cell-cycle: periodic sequence

• f cell divisions

crucial players: Cyclin, Cyclin-dependent kinase (Cdk) complicated regulatory network driven by an oscillator?

Mathematical modeling of cell cycle control networks ODE models Stochastic models Boolean models Network graphs Popular ODE models Mitotic oscillator (simple model of embryonic cell cycle) (Goldbeter, 1991) Yeast cell cycle (Chen, Tyson & Novak, 2004) Mammalian cell cycle (G´ erard & Goldbeter, 2012)

The mitotic oscillator has the following components Cyklin C active Cdk M inactive Cdk M+ active C-protease X inactive C-protease X+ C activates M+ ! M M activates X+ ! X X degrades C

• A. Goldbeter, PNAS (1991); more realistic larger models

contain subsystems similar to the Goldbeter model.

Dynamics is governed by Michaelis-Menten kinetics

cyklin C 0 aktive Cdk M 0 active Cyklin- protease X 0 ˙ C = vi vdX C Kd + C kdC ˙ M = V1 C Kc + C 1 M K1 + 1 M V2 M K2 + M ˙ X = V3M 1 X K3 + 1 X V4 X K4 + X

Michaelis constants Kj small Michaelis constants ε ⌧ 1 X ε + X = 8 < : ⇡ 1, X = O(1)

x 1+x, X = εx

⇡ 0, X = o(ε)

The mitotic oscillator has a periodic solution ˙ C = 1 4(1 X C) ˙ M = 6C 1 + 2C 1 M ε + 1 M 3 2 M ε + M ˙ X = M 1 X ε + 1 X 7 10 X ε + X parameter

kd = 0.25 vi = 0.25 Kc = 0.5 Kd = V1 = 3 V2 = 1.5 V3 = 1 V4 = 0.7

ε = 103

The periodic orbit lies in the cube [0, 1]3 ⇢ R3 partly very close to M = 0, X = 0, M = 1, X = 0

The periodic orbit lies in the cube [0, 1]3 ⇢ R3 partly very close to M = 0, X = 0, M = 1, X = 0 Theorem: ε ⌧ 1 ) exists periodic orbit Γε

Kosiuk + Sz. , J. Math.

• Bio. (2015)

singular perturbation?

The limit ε ! 0 contains a surprise For ε ! 0: mitotic oscillator limits on a discontinuous vector field sliding dynamics on sides and edges periodic orbit Γε converges to singular cycle Critical manifold: M = 0, M = 1 C = 0, C = 1 why? limε!0

X ε+X =

⇢ 1, X 6= 0 0, X = 0

Proof uses GSPT and blow-up Sliding on sides corresponds to slow motion on critical manifolds M = 0 and X = 0 M = 0 X = 0

Proof uses GSPT and blow-up Sliding on edge corresponds to slow motion on

• ne-dimensional critical manifold in blown-up edge.

• 4. Approach works / is useful for other small to

medium sized models negative feedback oscillator, Tyson et.al. (2003) mutual inhibition oscillator model of bipolar disorder, Goldbeter (2011) Repressilator, Elowitz & Leibler (2000) cell cycle models, Novak & Tyson (2001), Gerard & Goldbeter (2009, 2011) MAPK signaling, Kholodenko (2000)

Dynamics of Gerard & Goldbeter 39-variable model is “similar” to the Mitotic Oscillator GSPT analysis? start with their 5-variable Skeleton model

Slow-fast dynamics in models of MAPK-signaling Mitogen Activated Protein Kinase signaling pathway regulates cell growth, proliferation, cell death,... signal transduction from cell membrane to nucleus activation/deactivation: cascade of phosphorylation cycles

Mathematical models of the cascade based on mass action kinetics or Michaelis-Menten kinetics 9 variable model, Kholodenko (2000) 22 variable model, Huang & Ferrell (1996) Issues: signal amplifier (ultrasensitivity) bistability (biological switch)

• scillations (biological role mostly unclear)

feedbacks (some unknown or hypothetical) crosstalk to other pathways (some unknown or hypothetical) Existing mathematical work is based on numerical simulation.

Dynamics of Kholodenko’s MAPK-model is very similar to the Mitotic Oscillator MAPKpp solid line, MAPK dashed line Goal: analyze dynamics by GSPT based on small Michaelis-Menten constants

Interactions (crosstalk) with other pathways...

• 5. Many chemical, metabolic, and gene regulatory

networks show “nonstandard” or ‘hidden” slow-fast dynamics no global separation into slow and fast variables several scaling regimes with different limiting problems are needed singular or non-uniform dependence on several parameters, lack of smoothness need to combine different approaches: symbolic methods, tropical methods, reaction network theory, numerical simulation & bifurcation analysis, dynamical systems, GSPT and blow-up method,...

Our work suggests the following strategy identify fastest time-scale and corresponding scale of dependent variables, rescale

• ften the limiting problem has a (partially)

non-hyperbolic critical manifold use (repeated) blow-ups to desingularize identify relevant singular dynamics carry out perturbation analysis