Robustness of Complex Dynamics in Biology Workshop on Uncertain - - PDF document

• C. Breindl, D. Schittler, S. Waldherr, and F. Allgöwer

Institute for Systems Theory and Automatic Control University of Stuttgart, Germany

Robustness of Complex Dynamics in Biology

From qualitative to quantitative models

Workshop on Uncertain Dynamical Systems Udine, August 25, 2011

Uncertainty in models for biochemical networks

Biological knowledge

NF-κB I-κBα NF-κB I-κBα NF-κB I-κBα NF-κB I-κBα I-κBα mRNA nucleus cytosol IKK AAs

Experimental data ˙ x = F(x, p)

Types of uncertainty for biochemical networks

Structural: List of molecular species complete? Structural: Molecular interactions correct? Parametric: Parameter values correct?

Uncertainty

Robustness on different levels

˙ x = F(x, p)

Structural uncertainties Parametric uncertainties

A B C

? ? ? p ¯ x

Use assumption that nature favors robust solutions to "identify" likely interaction structure. How robust is the model’s behavior against perturbations in the parameters.

• F. Allgöwer, Robustness of complex dynamics in biology

Outline

1

Motivation

2

Steady state robustness of qualitative biological networks

3

Dynamical robustness of parametrized biological networks

4

Summary and conclusions

• F. Allgöwer, Robustness of complex dynamics in biology

Outline

1

Motivation

2

Steady state robustness of qualitative biological networks Qualitative modeling framework Definition and computation of robustness measure Application example: Differentiation of mesenchymal stem cells

3

Dynamical robustness of parametrized biological networks

4

Summary and conclusions

• F. Allgöwer, Robustness of complex dynamics in biology

Motivation

Typical situation when studying gene regulatory networks: List of species (genes, mRNAs, ...) given. Only little is known about the interactions between these species. Mostly only qualitative measurements are available. Two important questions:

1 Can a given interaction structure explain the observed behavior in principle? 2 Given several alternative model structures: Which one is biologically more

plausible?

• F. Allgöwer, Robustness of complex dynamics in biology

Important observed behavior: Multistability

Some example systems Life or death decision Differentiation

MSC tripotent

• steochondro-

progenitor pre-osteoblast

• steoblast

chondrocyte adipocyte fibroblast myoblast ...

Pattern formation

Prakash and Wurst, Cell.Mol.Life Sci. (2006)

Multistability is a widely observed phenomenon. How to model and analyze the system when only qualitative information about the steady state behavior is available?

• F. Allgöwer, Robustness of complex dynamics in biology

Modelling Framework: Network Structure

System equations: ˙ xi = −ki · xi + ϕi(x), i ∈ {1, . . . , n} ϕi(x) : arbitrary combination of sums and products of activation and inhibition functions Linear degradation A protein either activates or inhibits production of another protein ⇒ Monotonic activation- and inhibition functions

x1 x2

⇒

x1 act(x1) Activation function

˙ x2 = −k2 · x2 + act(x1)

x1 x2

⇒

x1 inh(x1) Inhibition function

˙ x2 = −k2 · x2 + inh(x1) Exact shapes of activation and inhibition functions not specified ⇒ structurally uncertain system

• F. Allgöwer, Robustness of complex dynamics in biology

Measurements

Steady state measurements Variability between systems Uncertain measurements ⇒ Representation of steady states as rectangular forward-invariant sets ⇒ Only distinction between “high” and “low” x1 x2 xmax

1

xmax

2

xhigh

2

xlow

2

xlow

1

xhigh

1

F2 F1 X

Mutual inhibition network: ˙ x1 = −k1 · x1 + inh1(x2) ˙ x2 = −k2 · x2 + inh2(x1)

Example for mutual inhibition network F1 = [0, xlow

1

] × [xhigh

2

, xmax

2

] = [0, 0.2] × [0.8, 1] F2 = [xhigh

1

, xmax

1

] × [0, xlow

2

] = [0.8, 1] × [0, 0.2]

• F. Allgöwer, Robustness of complex dynamics in biology

Outline Question 1: Can a given interaction structure explain the observed steady state behavior in principle?

• F. Allgöwer, Robustness of complex dynamics in biology

Specifications

Specification of (relative) concentrations and intervals 0 ≤ xlow

i

≤ xhigh

i

≤ xmax

i

for all i Ilow

xi

= [0, xlow

i

], Ihigh

xi

= [xhigh

i

, xmax

i

] Specification of m desired forward-invariant sets Fz = Ilow/high

x1

× . . . × Ilow/high

xn

z = 1, . . . , m Relates well to typical steady state measurements: Microarray data

• F. Allgöwer, Robustness of complex dynamics in biology

Specifications

Specification of (relative) concentrations and intervals 0 ≤ xlow

i

≤ xhigh

i

≤ xmax

i

for all i Ilow

xi

= [0, xlow

i

], Ihigh

xi

= [xhigh

i

, xmax

i

] Specification of m desired forward-invariant sets Fz = Ilow/high

x1

× . . . × Ilow/high

xn

z = 1, . . . , m Parametrization of activation/inhibition functions x inh(x) γlow γhigh N xlow xhigh xmax

• write ϕ T if ϕ lies in the tube T
• write ϕ T if ϕ violates tube T
• F. Allgöwer, Robustness of complex dynamics in biology

Conditions for forward-invariance

Forward-invariance of an invariant set F (Nagumo) Consider the system ˙ x = g(x). Assume that for each initial condition in a set X it admits a globally unique

• solution. Let F ⊆ X be a closed and convex set. Then F is

forward-invariant if and only if g(x) lies in the tangent cone to F in x for all x ∈ ∂F.

F

System satisfies assumptions. Every Fz is hyper-rectangular. Functions ϕ are monotonic. No nominal functions ϕ given x1 x2 xhigh

1

xmax

1

xlow

2

Fz

˙ x2 = −k2xlow

2

+ inh(x1)

• F. Allgöwer, Robustness of complex dynamics in biology

Conditions for forward-invariance

Forward-invariance of an invariant set F (Nagumo) Consider the system ˙ x = g(x). Assume that for each initial condition in a set X it admits a globally unique

• solution. Let F ⊆ X be a closed and convex set. Then F is

forward-invariant if and only if g(x) lies in the tangent cone to F in x for all x ∈ ∂F.

F

System satisfies assumptions. Every Fz is hyper-rectangular. Functions ϕ are monotonic. ⇒ Only vertices have to be considered. No nominal functions ϕ given x1 x2 xhigh

1

xmax

1

xlow

2

Fz

˙ x2 = −k2xlow

2

+ inh(x1)

• F. Allgöwer, Robustness of complex dynamics in biology

Conditions for forward-invariance

Forward-invariance of an invariant set F (Nagumo) Consider the system ˙ x = g(x). Assume that for each initial condition in a set X it admits a globally unique

• solution. Let F ⊆ X be a closed and convex set. Then F is

forward-invariant if and only if g(x) lies in the tangent cone to F in x for all x ∈ ∂F.

F

System satisfies assumptions. Every Fz is hyper-rectangular. Functions ϕ are monotonic. ⇒ Only vertices have to be considered. No nominal functions ϕ given ⇒ Reformulation using tubes as worst case approximation. x1 x2 xhigh

1

xmax

1

xlow

2

Fz

˙ x2 = −k2xlow

2

+ inh(x1)

γlow

• F. Allgöwer, Robustness of complex dynamics in biology

Conditions for forward-invariance

Invariance condition reformulated in terms of tubes

Given a set F = Il1

x1 × . . . × Iln xn with li ∈ {low, high} and given tubes T i,k that satisfy

the conditions ∀i ∈ {1, . . . n} : −ki · xi + γi,1 ◦ . . . ◦ γi,qi ≥ 0 (1) where xi = min

xi ∈Ili

xi

xi, and, with xj denoting the argument of ϕi,k, γi,k =

• if ϕi,k = νi,k ∧ 0 ∈ I

lj xj

min{γ : (x, γ) ∈ T i,k ∧ x ∈ I

lj xj }

• therwise

and ∀i ∈ {1, . . . n} : −ki · xi + γi,1 ◦ . . . ◦ γi,qi ≤ 0 (2) where xi = max

xi ∈Ili

xi

xi, and, with xj denoting the argument of ϕi,k, γi,k = max{γ : (x, γ) ∈ T i,k ∧ xj ∈ I

lj xj }.

If ∀i, k : ϕi,k T i,k, then the set F is forward-invariant.

• F. Allgöwer, Robustness of complex dynamics in biology

The structure validation problem

Consequence: Feasibility problem A given interaction structure can in principle explain the observed steady state behavior (forward-invariant sets) if and only if γ values (tubes) can be found such that Equations (1) and (2) hold. Often there are several model structures that can in principle explain the

• bserved behavior.

Can we find a measure for the biological plausibility of an interaction structure? Employ robustness considerations to develop this measure.

• F. Allgöwer, Robustness of complex dynamics in biology

Outline Question 2: Find a measure for the biological plausibility of a model structure based on robustness considerations.

• F. Allgöwer, Robustness of complex dynamics in biology

Robustness of multistable behavior

As before: Model should reproduce observed multistable behavior. Additionally: This behavior should be robust against perturbations in the interactions. Robustness of biological systems (Kitano, 2004) A system’s ability to maintain its function in the presence of perturbations. Here: System: Uncertain gene regulation network. Function: Existence of several forward-invariant sets at certain positions in the state space. Perturbations: Perturbations of the activation and inhibition functions, measured as P(ϕ, ϕp) = ∞ |ϕ(x) − ϕp(x)| dx. No nominal model given ⇒ Robustness measure for the interaction structure.

• F. Allgöwer, Robustness of complex dynamics in biology

The robustness measure

Mathematical formulation of robustness Minimal perturbation of ϕ to violate its tube T Rmin(ϕ, T) = min

ϕpT P(ϕ, ϕp)

ϕ that maximizes “robustness radius” (best centered ϕ) Rmax(T) = max

ϕT Rmin(ϕ, T)

x act(x) γlow γhigh γmax xlow xhigh xmax

Rmin

• F. Allgöwer, Robustness of complex dynamics in biology

The robustness measure

Mathematical formulation of robustness Minimal perturbation of ϕ to violate its tube T Rmin(ϕ, T) = min

ϕpT P(ϕ, ϕp)

ϕ that maximizes “robustness radius” (best centered ϕ) Rmax(T) = max

ϕT Rmin(ϕ, T)

x act(x) γlow γhigh γmax xlow xhigh xmax

Rmax Rmax

Robustness measure of the interaction structure R = max

T j min j

Rmax(T j) s.t. Invariance conditions (1) and (2) hold for every set Fz.

• F. Allgöwer, Robustness of complex dynamics in biology

Computation of the robustness measure

Rmax(T) can be computed analytically:

x act(x) γlow γhigh γmax xlow xhigh xmax

Rmax(T) = γlow · (γmax − γhigh) γlow + (γmax − γhigh) · (xhigh − xlow) Interpretation of R: Given optimal monotonous function ϕopt

j

with R = Ropt. Then each monotonous function ϕopt

j

can be perturbed by at least P(ϕopt

j

, ϕj) = Ropt and the sets Fz, z = 1, . . . , m, are still forward-invariant.

• F. Allgöwer, Robustness of complex dynamics in biology

Computation of the robustness measure

Rmax(T) can be computed analytically:

x act(x) γlow γhigh γmax xlow xhigh xmax

Rmax(T) = γlow · (γmax − γhigh) γlow + (γmax − γhigh) · (xhigh − xlow) If all monotonous functions are connected by “+”

• r all monotonous functions are connected by “·”

the problem of computing R is a convex optimization problem.

Details in C. Breindl et al., IJRNC, 2011

• F. Allgöwer, Robustness of complex dynamics in biology

Example: Differentiation of mesenchymal stem cells

Differentiation tree:

MSC tripotent

• steochondro-

progenitor pre-osteoblast

• steoblast

chondrocyte adipocyte fibroblast myoblast ...

• steochondro-

progenitor pre-osteoblast chondrocyte

Established idea: Hierarchical structure; binary differentiation tree Generalizes to other differentiation processes Modeling of one decision step

xC xO xP

Osteoblast Chondrocyte Progenitor Three cell types, three forward- invariant sets, three nodes. Interactions?

• F. Allgöwer, Robustness of complex dynamics in biology

Requirements on the model

Model structure p

• c

? ? ? Specification p

• c

progenitor high low low

• steoblast

low high low chondrocyte low low high Scan all possible 3 node networks

(using “·” for combined influences)

• F. Allgöwer, Robustness of complex dynamics in biology

Results

Building blocks

xi xi xi xi xi xi

every node xi, i ∈ {p, o, c}

• C. Breindl, D. Schittler et al., IFAC World Congress, 2011

Proof of principle: detailed analysis of one structure

xO xC xP zD zO zC

d dt xP = −kPxP + fP(xP, zD) d dt xO = −kOxO + fO(xO, xC, xP, zO) d dt xC = −kC xC + fC (xC, xO, xP, zC )

• D. Schittler et al., Chaos, 2010
• F. Allgöwer, Robustness of complex dynamics in biology

Results

Robustness versus number and sign of interconnections

1 2 3 4 5 6 7 8 9

#intercon. Robustness

−9 −7 −5 −3 −1 1 3 5 7 9

#act. - #inh. intercon.

p

• c

Maximally robust model

A higher number of interconnection tends to decrease robustness A highly robust network tends to have an excess of activating interconnections

• C. Breindl, D. Schittler et al., IFAC World Congress, 2011
• F. Allgöwer, Robustness of complex dynamics in biology

Results

Robustness versus maximal in-degree

• max. indeg.

Robustness

Maximal number of regulating factors of a gene limits robustness. Closer inspection of the results of the optimization problems: Nodes with most incoming interconnections form “bottlenecks”.

• C. Breindl, D. Schittler et al., IFAC World Congress, 2011
• F. Allgöwer, Robustness of complex dynamics in biology

Conclusions of Part I

Using an appropriate modeling framework for uncertain gene regulation networks, the presented method allows to answer the two questions posed:

1 We can check whether a given model structure can explain an observed

steady state behavior in principle.

2 We can evaluate how well a given model structure is suited to generate the

desired steady state behavior in terms of its robustness. Properties and advantages of the proposed robustness measure: The robustness measure R is a measure for the interaction structure itself: It indicates the robustness of the most robust realization within the given interaction structure. It can be used to compare and rank alternative model structures. It can be computed via a convex optimization problem for important classes of models. Efficient computational tool for robustness analysis in presence of only qualitative knowledge about the structure and the system behavior.

• F. Allgöwer, Robustness of complex dynamics in biology

Outline

1

Motivation

2

Steady state robustness of qualitative biological networks

3

Dynamical robustness of parametrized biological networks Robustness measure for dynamical behavior Feedback loop breaking Handelman representation and infeasibility certificates via Positivstellensatz Application example: Robustness of oscillations in the NF-κB pathway

4

Summary and conclusions

• F. Allgöwer, Robustness of complex dynamics in biology

Robustness of dynamical behavior

p ¯ x

Bistability

p ¯ x

Oscillations

Observations

Changes in the dynamical behavior occur at bifurcation points Complex dynamical behavior is coupled to unstable stationary points.

• F. Allgöwer, Robustness of complex dynamics in biology

The dynamical robustness radius

A robustness measure for dynamical behavior

The dynamical robustness radius is the smallest multiplicative parameter variation up to which no bifurcations occur.

(based on Ma & Iglesias 2002, extended to multi-parametric case) p1 p2 robust parameter set

Pr(ψ∗, p0) bifurcation surface

nominal p0 Pr(ψ, p0) =

• p ∈ Rm |

p0,i ψ

≤ pi ≤ ψp0,i, i = 1, . . . , m

• ⇒ dynamical robustness radius ψ∗
• F. Allgöwer, Robustness of complex dynamics in biology

Roles of feedback circuits

Positive feedback enables bistability (switching).

x1 x2 + +

Stimulus Response Time Protein activity

Negative feedback enables sustained oscillations.

x1 x2 x3 – – –

Time Protein activity

Observation

Feedback circuits enable complex dynamical behavior. Parameter values are relevant for behavior.

• F. Allgöwer, Robustness of complex dynamics in biology

Definition of feedback loop breaking

ODE model for biochemical network: ˙ x = Sv(x, p) = F(x, p)

Loop breaking definition

A loop breaking is a pair (f , h) such that F(x, p) = f (x, h(x), p) Open loop system: ˙ x = f (x, u, p) y = h(x) Closing the loop: u = h(x) ⇒ we recover the closed loop system

˙ x = f (x, u, p) y = h(x) y u

• F. Allgöwer, Robustness of complex dynamics in biology

Computing the feedback loop gain

Use feedback loop breaking to reformulate problem ˙ x = f (x, u, p) y = h(x) y u Linear approximation around stationary point (¯ x, p), F(¯ x, p) = 0, and Laplace transformation to frequency domain G(¯ x, p, s) y u

Re Im 1

G(¯ x, p, jω)

Result of reformulation

Transfer function G(¯ x, p, s) describes feedback loop gain Control theory (e.g. Nyquist criterion) can be used to check dynamical properties

• f the closed loop system
• F. Allgöwer, Robustness of complex dynamics in biology

Critical points: candidates for changes in the dynamical behavior

Definition: Critical point

The pair (¯ xc, pc), with F(¯ xc, pc) = 0, is called a critical point, if the Jacobian A(¯ xc, pc) = ∂F

∂x (¯

xc, pc) has an eigenvalue on the imaginary axis.

p ¯ x

Bistability Critical points

p ¯ x

Oscillations Critical points

• F. Allgöwer, Robustness of complex dynamics in biology

Characterization of critical points via the feedback loop gain

Theorem: Feedback loop gain and critical points

A critical point (¯ xc, pc) is characterized by the condition G(¯ xc, pc, jωc) = 1 Then, A(¯ xc, pc) has an eigenvalue at jωc.

Re Im

G(¯ xc, pc, jω)

1

ωc

Feedback loop gain G(¯ x, p, jω) is used to characterize candidate points for changes in the dynamical behavior.

• F. Allgöwer, Robustness of complex dynamics in biology

System of polynomial (in-)equalities

Do critical points exist?

G(x, p, jω) = 1 F(x, p) = 0 pi,min ≤ pi ≤ pi,max xj,min ≤ xj ≤ xj,max Specific problem: want to show that no solutions exist!

• F. Allgöwer, Robustness of complex dynamics in biology

Handelman representation

Consider a compact polytope K defined by the affine constraints Zj(x, p) ≥ 0, j = 1, . . . , M, Define Handelman monomials from the constraints Z(x, p): Hd(x, p) =

M

• j=1

Zj(x, p)dj , with d ∈ NM

0 the vectorial degree of the Handelman monomial Hd.

The polynomial Y : Rn+q → R is non-negative on K, if and only if Y can be represented as Y =

• d∈NM

cdHd, with non-negative coefficients cd.

• F. Allgöwer, Robustness of complex dynamics in biology

Positivstellensatz

Consider a system of polynomial (in-)equalities given by Yi(x, p) = 0, i = 1, . . . , N Zj(x, p) ≥ 0, j = 1, . . . , M. These equations do not have a solution, if and only if there exist Y in the ideal

• f Y1, . . . , YN and Z in the cone generated by Z1, . . . , ZM such that

Y + Z + 1 = 0. The ideal of Y1, . . . , YN is generated by N

i=1 TiYi.

The cone of Z1, . . . , ZM is generated by

d∈NM

0 cdHd.

• F. Allgöwer, Robustness of complex dynamics in biology

Positivstellensatz infeasibility certificates from linear programs

Taking again a system of polynomial (in-)equalities (PP) :

• Yi(x, p) = 0,

i = 1, . . . , N Zj(x, p) ≥ 0, j = 1, . . . , M, with affine Z.

By the Positivstellensatz:

(PP) does not have solutions ⇔       

N

• i=1

TiYi +

• d∈NM

cdHd + 1 = 0 cd ≥ 0 Infeasibility certificates for (PP) from a linear program in Ti and cd.

• F. Allgöwer, Robustness of complex dynamics in biology

Robustness certificates via feedback loop breaking

A robustness certificate: For a given parameter polytope P, try to guarantee that P is in the robust parameter set.

Robustness from the loop breaking characterization

G(¯ x, p, jω) = 1 F(¯ x, p) = 0 does not have a solution with p ∈

• P

⇔

• P ⊂ robust parameter set

Problem: Non-existence of solutions for system of polynomial equalities and affine inequalities Solution: New computational test with linear programming based on Positivstellensatz and Handelman refutation Computational test passed ⇒ Non-existence of solutions

Details in S. Waldherr and F. Allgöwer, Automatica, 2011

• F. Allgöwer, Robustness of complex dynamics in biology

Computation of the dynamical robustness radius

Dynamical robustness radius at nominal parameters p0 as

• ptimization problem:

ψ∗ = sup{ψ | Pr(ψ, p0) ⊂ robust parameter set} ⇒ combine robustness certificates with bisection on ψ

Pr(ψ∗, p0)

Pr(ψ3, p0)

1 Testing ψ1 . . . robust 2 Testing ψ2 . . . not robust 3 Testing ψ3 . . . robust

⇒ robustness radius ψ∗ ≥ ψ3

• F. Allgöwer, Robustness of complex dynamics in biology

Example: Application to NF-κB signaling

NF-κB I-κBα NF-κB I-κBα NF-κB I-κBα NF-κB I-κBα I-κBα mRNA nucleus cytosol IKK AAs

˙ x1 = ktl x2 − α(Ntot − x4)x1 KI + x1 − kI,inx1 + kI,outKNx3 KN + x4 ˙ x2 = ktx2

4 − γmx2

˙ x3 = kI,inx1 − kI,outKNx3 KN + x4 − kNI,outx3x4 KN + x4 ˙ x4 = kN,inKI (Ntot − x4) KI + x1 − kNI,outx3x4 KN + x4

400 800 t [min] nuclear NF-κB

• F. Allgöwer, Robustness of complex dynamics in biology

Robustness of NF-κB oscillations

Robustness radius computation

Loop breaking at I-κBα transcription Robustness certificates from linear program with 188489 parameters and 2282 equality constraints Robustness radius 2.078 < ψ∗ < 2.084

0.5 1 1.5 2 2.5 0.5 1 1.5

kt [1/(µM min)] α [1/min] Hopf bifurcation locus Region of certified robust instability p0

p0 ψ∗ ≤ p ≤ ψ∗p0 Oscillations are robust against significant parameter changes.

• F. Allgöwer, Robustness of complex dynamics in biology

Conclusions of Part II

Parametric uncertainty is an important issue when modeling biological systems. Feedback loop breaking useful to analyze dynamical properties. Analysis based on infeasibility certificates from polynomial programming, yields convex optimization problems. Efficient computational tool for robustness of small to medium scale biomolecular networks

Matlab toolbox for uncertainty analysis with polynomial programming

http://biosdp.sourceforge.net

• F. Allgöwer, Robustness of complex dynamics in biology

Summary and conclusions

Questions one wants to answer depend on the degree of uncertainty in the knowledge. Modeling and analysis tools must be able to deal with the type of uncertainty. Part I: A qualitative modeling and analysis approach to study how well a given interaction structure can reproduce an observed steady state behavior. Part II: A quantitative modeling and analysis approach to study the

• ccurrence of bifurcations in a given parameter region.

Both approaches guided by robustness considerations. Both approaches study truly nonlinear phenomena.

• F. Allgöwer, Robustness of complex dynamics in biology

Acknowledgment

Christian Breindl Daniella Schittler Steffen Waldherr

• F. Allgöwer, Robustness of complex dynamics in biology

References

1

• C. Breindl, D. Schittler, S. Waldherr, and F. Allgöwer

Structural requirements and discrimination of cell differentiation networks will be presented at 18th IFAC World Congress, 2011.

2

• C. Breindl, S. Waldherr, D.M. Wittmann, F.J. Theis, and F. Allgöwer

Steady-state robustness of qualitative gene regulation networks International Journal of Robust and Nonlinear Control, published online, 2011.

3

• D. Schittler, J. Hasenauer, S. Waldherr, and F. Allgöwer

Cell differentiation modeled via a coupled two-switch regulatory network Chaos (20), art.no. 045121, 2010.

4

• S. Waldherr and F. Allgöwer

Robust stability and instability of biochemical networks with parametric uncertainty Automatica 47(6): 1139-1146, 2011

• F. Allgöwer, Robustness of complex dynamics in biology