From RNA-seq Time Series Data to Models of Regulatory Networks - - PowerPoint PPT Presentation

From RNA-seq Time Series Data to Models of Regulatory Networks

Konstantin Mischaikow

• Dept. of Mathematics, Rutgers

mischaik@math.rutgers.edu

Brown, Dec 2017

A Model Problem

Estimated number of malaria cases in 2010: between 219 and 550 million Estimated number of deaths due to malaria in 2010: 600,000 to 1,240,000 Malaria may have killed half of all the people that ever lived. And more people are now infected than at any point in history. There are up to half a billion cases every year, and about 2 million deaths - half of those are children in sub-Saharan Africa.

• J. Whitfield, Nature, 2002

Resistance is now common against all classes of antimalarial drugs apart from artemisinins. … Malaria strains found on the Cambodia– Thailand border are resistant to combination therapies that include artemisinins, and may therefore be untreatable. World Health Organization Malaria is of great public health concern, and seems likely to be the vector-borne disease most sensitive to long-term climate change. World Health Organization

Malaria

Malaria: P. falciparum

48 hour cycle 1-2 minutes Task: Understand the regulation on the genetic/biomolecular level with the goal of affecting the dynamics with drugs. A proposed network A differential equation dx

dt = f(x, λ) is proba-

bly a reasonable model for the dynamics, but I do not have an analytic description of f or estimates of the parameters λ. Malaria is

• Sequenced
• Poorly annotated

All genes (5409)

1.5$ 0.0$ &1.5$ Standard$devia0ons$from$ mean$expression$(z&score)$

High Low 0$ 10$ 20$ 30$ 40$ 0me$in#vitro#(hours)$ 50$ 60$ periodic$genes$(43)$

10 20 30 40 50 60

Walter Reed Army Inst. Research Duke University

A Philosophical Intfrlude

The Lac Operon

Ozbudak et al. Nature 2004

Network Model 1 τy ˙ y = α RT RT + R(x) − y 1 τx ˙ x = βy − x R(x) = RT 1 + ⇣

x x0

⌘n ODE Model Data ODES are great modeling tools, but should be handled with care. parameter values

α = 84.4 1 + (G/8.1)1.2 + 16.1 β = . . .

Classical Qualita7ve Representa7on

• f Dynamics

Dynamic Signature (Morse Graph)

Not Precise Accurate Rigorous Precise Not Accurate Not Rigorous What does it mean to solve an ODE? Conley-Morse Chain Complex

model “truth” parameter

Combinatprial Dynamics

p1 p0 p2 p3

Vertices: States Edges: Dynamics

Don’t know exact current state, so don’t know exact next state

Simple decomposition

• f Dynamics:

Recurrent Nonrecurrent (gradient-like)

Linear time Algorithm!

Morse Graph

• f state transition graph

State Transition Graph

F : X − → → X

Poset

p1 p0 p2 p3

P O S E T

Morse Graph

• f F : X −

→ → X Join Irreducible J∨(A)

Birkhoff’s Theorem implies that the Morse graph and the lattice

• f Attractors are equivalent.

What is observable?

A ⊂ X is an attractor if F(A) = A

p1 p0 p1, p0 p2, p1, p0 p3, p2, p1, p0 Lower Sets O(M)

Computable Observable

∅ Lattice of Attractors

• f F : X −

→ → X ∨ = ∪ ∧ = maximal attractor in ∩

Topology

(differentjal equatjons are not defined on discretf sets)

Let X be a compact metric space. phase space Let R(X) denote the lattice of regular closed subsets of X. Infinite unbounded lattice Level of measurement Applicable scale for model Let L be a finite bounded sublattice

• f R(X).

Topology Dynamics Use Birkhoff to define poset (P := J∨(A), <) Remark: I have purposefully ignored the relation between L and F : X − → → X G(L) denoted atoms of L “smallest” elements of L For each p ∈ P define a Morse tile M(p) := cl(A \ pred(A)) Declare a bounded sublattice A ⊂ L to be a lattice of forward invariant regions (attracting blocks).

Example

Morse tiles M(p) Let F 0(x) = −f(x).

• 4

4

Atoms of lattice: G(L) = {[n, n + 1] | n = −4, . . . , 3} Phase space: X = [−4, 4] ⊂ R P

1 2 3

Birkhoff

How does this relate to a differential equation dx

dt = f(x)?

• 4

4

F F A Lattice of attracting blocks: A = {[−3, −1], [1, 3], [−3, −1] ∪ [1, 3], [−4, 4]} Attracting blocks are regions of phase space that are forward invariant with time. Remark: This leads to a homology theory F

Switching Systfms

(an example of how tp use tiese ideas)

Choosing L and F : X − → → X

Biological Model

Assume xi decays.

dxi dt = −γixi dxi dt = −γixi + Λi(x) dxi dt = −γixi + Λi(xj)

How do I want to interpret this information? What differential equation do I want to use? Proposed model:

dx2 dt x1 θ2,1 u2,1 l2,1

x1 represses the production of x2.

1 2

x1 activates the production of x2.

1 2

Parameters 1/node 3/edge For x1 < θ2,1 we ask about sign (−γ2x2 + u2,1). For x1 > θ2,1 we ask about sign (−γ2x2 + l2,1). xi denotes amount of species i.

−

j,i(xi) =

( uj,i if xi < ✓j,i `j,i if xi > ✓j,i

Focus on sign of −γixi + σ−

i,j(xj)

• −γixi + σ+

i,j(xj)

x1 x2

Phase space: X = (0, ∞)2 Parameter space is a subset of (0, ∞)8

θ2,1 θ1,2

Fix z a regular parameter value. z is a regular parameter value if 0 < i 0 < `i,j < ui,j, 0 < ✓i,k 6= ✓j,k, and 0 6= i✓j,i + Λi(x)

Example

1 2

−γ2θ1,2 + σ+

2,1(x1)

−γ1θ2,1 + σ−

1,2(x2)

We care about sign of: If x1 < θ2,1 and −γ2θ1,2 + σ+

2,1(x1) > 0

If x1 < θ2,1 and −γ2θ1,2 + σ+

2,1(x1) < 0

θ2,1 θ1,2 x1 x2

Need to Construct State Transition Graph Fz : X − → → X Fix z a regular parameter value. Vertices X corresponds to all rectangular domains and co-dimension 1 faces defined by thresholds θ. Faces pointing in map to their domain. Domains map to their faces pointing

• ut.

Edges If no outpointing faces domain maps to itself.

θ2,1 θ1,2 x1 x2

Morse Graph FP{0,1} Fix z a regular parameter value.

Example

1 2 θ2,1 θ1,2 x1 x2

Assume: l1,2 < γ1θ2,1 < u1,2 γ2θ1,2 < l2,1 < u2,1 l1,2 < γ1θ2,1 < u1,2 l2,1 < γ2θ1,2 < u2,1 Check signs of −γiθj,i + σ±

i,j(xj)

Constructing state transition graph Fz : X − → → X (M1, M2, m1, m2) Dynamics orders maxima and minima FC

DSGRN Database Parameter graph provides explicit partition of entire 8-D parameter space. Observe that we can query this database for local or global dynamics. Input: Regulatory Network

1 2

Output: DSGRN database

(1) FP(1,1) γ1θ2,1 < l1,2 < u1,2 γ2θ1,2 < l2,1 < u2,1 (4) FP(1,1) γ1θ2,1 < l1,2 < u1,2 l2,1 < γ2θ1,2 < u2,1 (7) FP(1,0) γ1θ2,1 < l1,2 < u1,2 l2,1 < u2,1 < γ2θ1,2 (2) FP(0,1) l1,2 < γ1θ2,1 < u1,2 γ2θ1,2 < l2,1 < u2,1 (5) FC l1,2 < γ1θ2,1 < u1,2 l2,1 < γ2θ1,2 < u2,1 (8) FP(1,0) l1,2 < γ1θ2,1 < u1,2 l2,1 < u2,1 < γ2θ1,2 (3) FP(0,1) l1,2 < u1,2 < γ1θ2,1 γ2θ1,2 < l2,1 < u2,1 (6) FP(0,0) l1,2 < u1,2 < γ1θ2,1 l2,1 < γ2θ1,2 < u2,1 (9) FP(0,0) l1,2 < u1,2 < γ1θ2,1 l2,1 < u2,1 < γ2θ1,2

Back tp Malaria

Remark: there are a variety of statistical methods for generating possible regulatory networks from this type of time series data.

In vitro data (WRAIR+ Haase)

High Low 1.5

• 1.5

Standard deviations from mean expression (z-score) 10 20 30 40 50 60 Time in vitro (hours) Putative TF genes (456)

Remarks about dynamics:

• 1. Gene expression is

cyclic in nature.

• 2. We know relative

times of expression

• f genes.

Assumption: Expression of important functions must be robust to perturbations.

Experimental 7me series for associated genes

Simple Test

Cyclic feedback system: well understood using classical dynamical systems techniques. Under the assump7on of monotone switches if parameter values are chosen such that there exists a stable periodic orbit, then the maxima in the network must occur in the order: (188,93,184, 395) (green, blue, cyan,red)

Conclusion: This network does not generate observed dynamics

No mathema7cal theory Time series for associated genes Computa7on 7me on laptop approximately 1 second. SQL Query: A stable cycle involving

• scilla7ons in all genes

DSGRN computa7on produces a parameter graph with approximately 45,000 nodes. 96 parameter graph nodes with Morse graph that has a minimal node consis7ng of a Full Cycle (FC).

Have developed polynomial 7me algorithm that take paths in state transi7on graph and iden7fies sequences of possible maxima and minima. Tested all max-min sequences from state transi7on graphs from all 96 parameter graph nodes against 17,280 experimental pa`erns. No Match

Conclusion: This network does not generate observed dynamics

DSGRN Analysis (II): Max-Min Matching

M m M m M m M M m M m m

DSGRN strategy

• I. Start with a proposed gene regulatory network.

Perform DSGRN computation to identify parameter node for which minimal Morse node is FC Reject parameter node if max-min sequences of FC are not linear extensions of poset. Extract from experimental data the poset indicating possible max-min orderings. Compute fraction of parameter nodes that match experimental data. IF fraction is small, THEN reject regulatory network. Assumption: Expression of important functions must be robust to perturbations.

DSGRN strategy

• II. Start with an acceptable gene regulatory network.

Add/Remove edge(s) Create new regulatory network via random perturbations: Add/Remove node

• I. Start with the newly proposed gene

regulatory network.

Current favorite network

90% of parameter nodes* induce minimal FC node 80% of parameter nodes* induce minimal FC which agrees with experimentally determined max-min

• rdering.

Thank-you for your Attention

Homology + Database Software chomp.rutgers.edu

Rutgers

• S. Harker

MSU

• T. Gedeon
• B. Cummings

FAU

• W. Kalies

VU Amsterdam

• R. Vandervorst

Example (continued)

• 4

4

X = [−4, 4] ⊂ R

1

The homology Conley index of M(p) is CH∗(p) := H∗(A, pred(A); k) k a field

2 3

(P, <) (k, 0, . . .) (k, 0, . . .)   1 1   A Theorem: (R. Franzosa) There exists a strictly upper triangular (with respect to <) boundary operator ∆: M

p∈P

CH∗(p) → M

p∈P

CH∗(p) such that the induced homology is isomorphic to H∗(X). (0, k, 0, . . .) Conley index can be used to guarantee existence of equilibria, periodic

• rbits, heteroclinic and homoclinic orbits, and chaotic dynamics.

F

F