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

Malaria 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 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 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: P. falciparum 48 hour cycle Task: Understand the regulation on the genetic/biomolecular level with the goal of affecting the dynamics 1-2 minutes with drugs. All genes (5409) A differential equation dx dt = f ( x, λ ) is proba- High 1.5$ bly a reasonable model for the dynamics, but mean$expression$(z&score)$ Standard$devia0ons$from$ periodic$genes$(43)$ I do not have an analytic description of f or 0.0$ estimates of the parameters λ . &1.5$ Low Malaria is 0$ 10$ 20$ 30$ 40$ 50$ 60$ 0me$ in#vitro# (hours)$ • Sequenced • Poorly annotated 0 10 20 30 40 50 60 Walter Reed Army Inst. Research A proposed network Duke University

A Philosophical In tf rlude

The Lac Operon Ozbudak et al. Nature 2004 ODES are great modeling tools, but should be handled with care. Network Model 1 R T y = α ˙ R T + R ( x ) − y τ y 1 Data x = β y − x ˙ τ x 84 . 4 R T α = 1 + ( G/ 8 . 1) 1 . 2 + 16 . 1 R ( x ) = ⌘ n ⇣ x β = . . . 1 + x 0 parameter values ODE Model

What does it mean to solve an ODE? “truth” model Precise Not Accurate Not Rigorous Classical Qualita7ve parameter Representa7on of Dynamics Conley-Morse Chain Complex Not Precise Accurate Rigorous Dynamic Signature (Morse Graph)

Combina tp rial Dynamics

State Transition Graph F : X − → → X p 3 Linear time Algorithm! Simple decomposition of Dynamics: p 2 Recurrent Nonrecurrent (gradient-like) p 0 p 1 Morse Graph Vertices: States of state transition graph Edges: Dynamics Don’t know exact current state, Poset so don’t know exact next state

What is observable? A ⊂ X is an attractor if F ( A ) = A Birkhoff’s Theorem implies that the Morse graph and the lattice p 3 of Attractors are equivalent. Observable Computable P p 3 , p 2 , p 1 , p 0 O Lower Sets O ( M ) S p 2 , p 1 , p 0 p 2 E Join Irreducible T J ∨ ( A ) p 1 , p 0 p 1 p 0 p 0 p 1 ∅ Lattice of Attractors Morse Graph ∨ = ∪ of F : X − → → X ∧ = maximal attractor in ∩ of F : X − → → X

Topology ( di ff eren tj al equa tj ons are not de fi ned on discre tf sets )

Let X be a compact metric space. phase space Topology Let R ( X ) denote the lattice of Infinite unbounded regular closed subsets of X . lattice Let L be a finite bounded sublattice Level of measurement of R ( X ) . Applicable scale for model G ( L ) denoted atoms of L “smallest” elements of L Declare a bounded sublattice A ⊂ L to be a lattice of forward invariant regions (attracting blocks). Dynamics Use Birkhoff to define poset ( P := J ∨ ( A ) , < ) For each p ∈ P define a Morse tile M ( p ) := cl( A \ pred ( A )) Remark: I have purposefully ignored the relation between L and F : X − → → X

Example Phase space: X = [ − 4 , 4] ⊂ R Atoms of lattice: G ( L ) = { [ n, n + 1] | n = − 4 , . . . , 3 } Lattice of attracting blocks: A = { [ − 3 , − 1] , [1 , 3] , [ − 3 , − 1] ∪ [1 , 3] , [ − 4 , 4] } Morse tiles M ( p ) 3 Birkhoff P A -4 0 4 1 2 F How does this relate to a differential equation dx dt = f ( x ) ? F Let F 0 ( x ) = − f ( x ) . F Attracting blocks are regions of phase space that are forward invariant with time. -4 0 4 Remark: This leads to a homology theory

Switching Sys tf ms ( an example of how tp use ti ese ideas ) Choosing L and F : X − → → X

Biological Model How do I want to interpret this information? What differential equation do I want to use? x i denotes amount of species i . dx i Assume x i decays. Parameters dt = − γ i x i 1/node 3/edge dx i dx i Proposed model: dt = − γ i x i + Λ i ( x ) dt = − γ i x i + Λ i ( x j ) dx 2 1 2 dt 1 2 u 2 , 1 x 1 represses the x 1 activates the production of x 2 . production of x 2 . l 2 , 1 x 1 θ 2 , 1 For x 1 < θ 2 , 1 we ask about sign ( − γ 2 x 2 + u 2 , 1 ) . ( if x i < ✓ j,i u j,i j,i ( x i ) = � − if x i > ✓ j,i ` j,i For x 1 > θ 2 , 1 we ask about sign ( − γ 2 x 2 + l 2 , 1 ) . � − γ i x i + σ + � Focus on sign of − γ i x i + σ − i,j ( x j ) i,j ( x j )

− γ 1 θ 2 , 1 + σ − 1 , 2 ( x 2 ) Example We care about sign of: 1 2 − γ 2 θ 1 , 2 + σ + 2 , 1 ( x 1 ) Parameter space is a subset of (0 , ∞ ) 8 Fix z a regular parameter value. x 2 If x 1 < θ 2 , 1 and − γ 2 θ 1 , 2 + σ + 2 , 1 ( x 1 ) > 0 If x 1 < θ 2 , 1 and − γ 2 θ 1 , 2 + σ + 2 , 1 ( x 1 ) < 0 θ 1 , 2 z is a regular parameter value if 0 < � i 0 < ` i,j < u i,j , x 1 θ 2 , 1 0 < ✓ i,k 6 = ✓ j,k , and Phase space: X = (0 , ∞ ) 2 0 6 = � � i ✓ j,i + Λ i ( x )

Fix z a regular parameter value. Need to Construct State Transition Graph F z : X − → → X Vertices X corresponds to all rectangular x 2 domains and co-dimension 1 faces defined by thresholds θ . Edges θ 1 , 2 Faces pointing in map to their domain. Domains map to their faces pointing out. x 1 If no outpointing faces domain maps θ 2 , 1 to itself.

Example Constructing state transition graph F z : X − → → X 1 2 Check signs of − γ i θ j,i + σ ± i,j ( x j ) Fix z a regular parameter value. Assume: l 1 , 2 < γ 1 θ 2 , 1 < u 1 , 2 l 1 , 2 < γ 1 θ 2 , 1 < u 1 , 2 γ 2 θ 1 , 2 < l 2 , 1 < u 2 , 1 l 2 , 1 < γ 2 θ 1 , 2 < u 2 , 1 x 2 x 2 θ 1 , 2 θ 1 , 2 x 1 x 1 θ 2 , 1 θ 2 , 1 Dynamics orders Morse FP{0,1} FC maxima and minima Graph ( M 1 , M 2 , m 1 , m 2 )

DSGRN Database Output: Input: 1 2 DSGRN database Regulatory Network FP(1,0) FP(1,0) FP(0,0) γ 1 θ 2 , 1 < l 1 , 2 < u 1 , 2 l 1 , 2 < γ 1 θ 2 , 1 < u 1 , 2 l 1 , 2 < u 1 , 2 < γ 1 θ 2 , 1 (7) (8) (9) l 2 , 1 < u 2 , 1 < γ 2 θ 1 , 2 l 2 , 1 < u 2 , 1 < γ 2 θ 1 , 2 l 2 , 1 < u 2 , 1 < γ 2 θ 1 , 2 FP(1,1) FC FP(0,0) γ 1 θ 2 , 1 < l 1 , 2 < u 1 , 2 l 1 , 2 < γ 1 θ 2 , 1 < u 1 , 2 l 1 , 2 < u 1 , 2 < γ 1 θ 2 , 1 (4) (5) (6) l 2 , 1 < γ 2 θ 1 , 2 < u 2 , 1 l 2 , 1 < γ 2 θ 1 , 2 < u 2 , 1 l 2 , 1 < γ 2 θ 1 , 2 < u 2 , 1 FP(1,1) FP(0,1) FP(0,1) γ 1 θ 2 , 1 < l 1 , 2 < u 1 , 2 l 1 , 2 < γ 1 θ 2 , 1 < u 1 , 2 l 1 , 2 < u 1 , 2 < γ 1 θ 2 , 1 (1) (2) (3) γ 2 θ 1 , 2 < l 2 , 1 < u 2 , 1 γ 2 θ 1 , 2 < l 2 , 1 < u 2 , 1 γ 2 θ 1 , 2 < l 2 , 1 < u 2 , 1 Parameter graph provides explicit partition of entire 8 -D parameter space. Observe that we can query this database for local or global dynamics.

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) Remarks about dynamics: High 1. Gene expression is 1.5 cyclic in nature. Putative TF genes (456) Standard deviations from mean expression (z-score) 2. We know relative times of expression of genes. 0 -1.5 Low Assumption: Expression of 10 20 30 40 50 0 60 important functions must be robust to perturbations. Time in vitro (hours)

Simple Test Cyclic feedback system: Experimental 7me series for well understood using associated genes 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 DSGRN computa7on produces a SQL Query: A stable cycle involving parameter graph with oscilla7ons in all genes approximately 45,000 nodes. 96 parameter graph nodes with Computa7on 7me on laptop Morse graph that has a minimal approximately 1 second. node consis7ng of a Full Cycle (FC).

DSGRN Analysis (II): Max-Min Matching M m m M m M M M m M m m 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