ODEs Based Signaling Pathways Dynamics P.S. Thiagarajan School of - - PowerPoint PPT Presentation

odes based signaling pathways
SMART_READER_LITE
LIVE PREVIEW

ODEs Based Signaling Pathways Dynamics P.S. Thiagarajan School of - - PowerPoint PPT Presentation

Probabilistic Approximations of ODEs Based Signaling Pathways Dynamics P.S. Thiagarajan School of Computing, National University of Singapore Biopathways Biopathways: Metabolic Pathways Signaling Pathways Gene Regulatory


slide-1
SLIDE 1

Probabilistic Approximations of ODEs Based Signaling Pathways Dynamics

P.S. Thiagarajan

School of Computing, National University of Singapore

slide-2
SLIDE 2

Biopathways

 Biopathways:

 Metabolic Pathways  Signaling Pathways  Gene Regulatory Networks

slide-3
SLIDE 3

Signaling Pathways

  • Chemical reactions in

response to external signals (ligands)

  • Signals pass into the nucleus

through a series of protein modifications „Data transfer‟ mechanism of the cell

slide-4
SLIDE 4

A Common Modeling Approach

  • View a pathway as a network of bio-chemical

reactions

  • Model the network as a system of ODEs
  • One for each molecular species
  • Reaction kinetics: Mass action law, Michelis-Menten,

Hill, etc.

  • Study the ODE system dynamics.
slide-5
SLIDE 5

Assume mass law.

The ODEs model

฀  dS dt  k1  S  E  k2  ES dE dt  k1  S  E  (k2  k3) ES dES dt  k1  S  E  (k2  k3) ES dP dt  k3  ES

ES E S  P E 

฀  k1

฀  k3

฀  k2

slide-6
SLIDE 6

 Alternative approach:

 Keep track of exact number of molecules of

each type. Simulate the dynamics by executing one reaction at a time stochastically (CTMCs)

Stochastic simulations (Gillespie‟s algorithm)  Kappa , BioNetGen, PRISM, Bio-Pepa, ..

slide-7
SLIDE 7

ODEs: Major Hurdles

  • Many unknown rate constants.
  • Must be estimated using limited data:
  • Low precision, population-based, noisy
slide-8
SLIDE 8

Major Hurdles

  • High dimensional non-linear system
  • no closed-form solutions
  • must resort to numerical simulations
  • point values of initial states/data will not be available
  • a large number of numerical simulations needed for

answering each analysis question

slide-9
SLIDE 9

“Polling” based approximation

  • Start with an ODEs system.
  • Discretize the time and value domains.
  • Assume a (uniform) distribution of initial states
  • Generate a “sufficiently” large number of

trajectories by

  • Sampling the initial states and numerical simulations.
slide-10
SLIDE 10

The “exit poll” Idea

  • Encode this collection of discretized trajectories

as a dynamic Bayesian network.

  • ODEs  DBN
  • Pay the one-time cost of constructing the DBN

approximation.

  • Do analysis using Bayesian inferencing on the

DBN.

slide-11
SLIDE 11

Time Discretization

  • Observe the system only at a finite number of

time points.

x(t) t0 t1 t2 tmax ... ... ... ... x(t) = t3 + 4t + 2 2 x(0) = 2

slide-12
SLIDE 12

Value Discretization

  • Observe only with bounded precision

x(t) t0 t1 t2 tmax ... ... ... ... A B C D E

slide-13
SLIDE 13

x(t) t0 t1 t2 tmax ... ... ... ... A B C D E

  • A trajectory is recorded as a finite sequence of

discrete values.

Symbolic trajectories

(C,0) (D,1) (D,2) (E,3) (D,4) (C,5) (B,6)

slide-14
SLIDE 14

Collection of Trajectories

  • Assume a prior distribution of the initial states.
  • Uncountably many trajectories. Represented as a set of

(timed) finite sequences.

t x(t) t0 t1 t2 tmax ... ... ... ... A B C D E ... ...

(C,0) (D,1) (D,2) (E,3) (D,4) (C,5) (B,6) (D,3)

slide-15
SLIDE 15

Piecing trajectories together..

  • In fact, a probabilistic transition system.
  • Pr( (D, 2) (E, 3) ) is

the “fraction” of the trajectories residing in D at t = 2 that land in E at t = 3.

t t0 t1 t2 tmax ... ... ... ... A B C D E ... ...

(C,0) (D,1) (D,2) (E,3) (D,4) (C,5) (B,6) (D,3) 0.8 0.2

slide-16
SLIDE 16

The Justification

  • The value space of the variables is assumed to be a

compact subset C of

  • In Z’ = F(Z), F is assumed to be continuously

differentiable in C.

  • Mass-law, Michaelis-Menton,…
  • Then the solution t : C  C (for each t) exists, is

unique, a bijection, continuous and hence measurable.

  • But the transition probabilities can’t be computed.
slide-17
SLIDE 17

(s, i) – States; (s, i)  (s‟, i+1) -- Transitions Sample, say, 1000 times the initial states. Through numerical simulation, generate 1000 trajectories. Pr((s, i)  (s‟ i+1)) is the fraction of the trajectories that are in s at ti which land in s‟ at t i+1. t t0 t1 t2 tmax ... ... ... ... A B C D E

A computational approximation

1000 800

... ...

(C,0) (D,1) (D,2) (E,3) (D,4) (C,5) (B,6) (D,3) 0.8 0.2

slide-18
SLIDE 18

Infeasible Size!

  • But the transition system will be huge.
  • O(T . kn)
  • k  2 and n ( 50-100).
slide-19
SLIDE 19

Compact Representation

  • Exploit the network structure (additional

independence assumptions) to construct a DBN instead.

  • The DBN is a factored form of the probabilistic

transition system.

slide-20
SLIDE 20

Assume mass law.

The DBN representation

฀  dS dt  k1  S  E  k2  ES dE dt  k1  S  E  (k2  k3) ES dES dt  k1  S  E  (k2  k3) ES dP dt  k3  ES

ES E S  P E 

฀  k1

฀  k3

฀  k2

slide-21
SLIDE 21

฀  dS dt  k1  S  E  k2  ES dE dt  k1  S  E  (k2  k3) ES dES dt  k1  S  E  (k2  k3) ES dP dt  k3  ES

ES E S  P E 

฀  k1

฀  k3

฀  k2

S E ES P Dependency diagram

slide-22
SLIDE 22

Dependency diagram

฀  dS dt  k1  S  E  k2  ES dE dt  k1  S  E  (k2  k3) ES dES dt  k1  S  E  (k2  k3) ES dP dt  k3  ES

ES E S  P E 

฀  k1

฀  k3

฀  k2

S E ES P

slide-23
SLIDE 23

The DBN Representation

S0 E0 ES0 P0 S1 E1 ES1 P1 S2 E2 ES2 P2 S3 E3 ES3 P3

... ... ... ... ... ... ... ... ฀  dS dt  k1  S  E  k2  ES dE dt  k1  S  E  (k2  k3) ES dES dt  k1  S  E  (k2  k3) ES dP dt  k3  ES

ES E S  P E 

฀  k1

฀  k3

฀  k2

slide-24
SLIDE 24
  • Each node has a CPT

associated with it.

  • This specifies the local

(probabilistic) dynamics.

S0 E0 ES0 P0 P1 E2 ES2 P2 S3 E3 ES3 P3

... ... ... ... ... ... ... ... A B C

S1 E1 ES1 S2 P(S2=C|S1=B,E1=C,ES1=B)= 0.2 P(S2=C|S1=B,E1=C,ES1=C)= 0.1 P(S2=A|S1=A,E1=A,ES1=C)= 0.05

. . .

slide-25
SLIDE 25
  • Fill up the entries in

the CPTs by sampling, simulations and counting

S0 E0 ES0 P0 P1 E2 ES2 P2 S3 E3 ES3 P3

... ... ... ... ... ... ... ... A B C

S1 E1 ES1 S2

P(S2=C|S1=B,E1=C,ES1=B)= 0.2 P(S2=C|S1=B,E1=C,ES1=C)= 0.1 P(S2=A|S1=A,E1=A,ES1=C)= 0.05 . . .

slide-26
SLIDE 26

Computational Approximation

  • Fill up the entries in

the CPTs by sampling, simulations and counting

S0 E0 ES0 P0 P1 E2 ES2 P2 S3 E3 ES3 P3

... ... ... ... ... ... ... ... A B C

S1 E1 ES1 S2

500 100 1000

slide-27
SLIDE 27

The Technique

  • Fill up the entries in

the CPTs by sampling, simulations and counting

S0 E0 ES0 P0 P1 E2 ES2 P2 S3 E3 ES3 P3

... ... ... ... ... ... ... ... A B C

S1 E1 ES1 S2

500 100 P(S2=C|S1=B,E1=C,ES1=B)= 100/500= 0.2

slide-28
SLIDE 28

The Technique

S0 E0 ES0 P0 P1 E2 ES2 P2 S3 E3 ES3 P3

... ... ... ... ... ... ... ...

S1 E1 ES1 S2

The size of the DBN is: O(T . n . kd) d will be usually much smaller than n.

slide-29
SLIDE 29

Unknown rate constants

฀  dS dt  0.1 S  E  0.2 ES dE dt  0.1 S  E  (0.2  k3) ES dES dt  0.1 S  E  (0.2  k3) ES dP dt  k3  ES

ES E S  P E 

฀  k1  0.1

฀  k3

฀  k2  0.2

S0 E0 ES0 P0 S1 E1 ES1 P1 S2 E2 ES2 P2 S3 E3 ES3 P3

... ... ... ... ... ... ... ...

k3

dt dk3 = 0

slide-30
SLIDE 30

Unknown rate constants

During the numerical generation of a trajectory, the value

  • f k3 does not change

after sampling.

S0 E0 ES0 P0 S1 E1 ES1 P1 S2 E2 ES2 P2 S3 E3 ES3 P3

... ... ... ... ... ... ... ...

k3 k3 k3 k3

... ... P(k3=A|k3=A)= 1

1 2 3 1

slide-31
SLIDE 31

Unknown rate constants

During the numerical generation of a trajectory, the value

  • f k3 does not change

after sampling.

S0 E0 ES0 P0 S1 E1 ES1 P1 S2 E2 ES2 P2 S3 E3 ES3 P3

... ... ... ... ... ... ... ...

k3 k3 k2 k3

... ...

1 2 3

P(k3=A|k3=A)= 1 P(ES2=A|S1=C,E1=B,ES1=A,k3=A)= 0.4

1 1

slide-32
SLIDE 32

Unknown rate constants

Sample uniformly across all the Intervals.

S0 E0 ES0 P0 S1 E1 ES1 P1 S2 E2 ES2 P2 S3 E3 ES3 P3

... ... ... ... ... ... ... ...

k3 k3 k3 k3

... ...

1 2 3

slide-33
SLIDE 33

DBN based Analysis

  • Use Bayesian inferencing to do parameter estimation,

sensitivity analysis, probabilistic model checking …

  • Exact inferencing is not feasible for large models.
  • We do approximate inferencing.
  • Factored Frontier algorithm.
slide-34
SLIDE 34

Parameter Estimation

S0 E0 ES0 P0 S1 E1 ES1 P1 S2 E2 ES2 P2 S3 E3 ES3 P3

... ... ... ... ... ... ... ...

k3 k3 k3 k3

  • 1. For each choice of (interval)

values for unknown parameters, run FF, compare with experimental data and assign a score using FF.

  • 2. Return parameter estimates as

maximal likelihoods.

  • 3. FF can be then used on the

calibrated model to do sensitivity analysis, probabilistic verification etc. ... ...

1 2 3

slide-35
SLIDE 35

DBN based Analysis

  • Our experiments with signaling pathways

models (taken from the BioModels data base) show:

  • The one-time cost of constructing the DBN can be

easily amortized by using it to do parameter estimation and sensitivity analysis.

  • Good compromise between efficiency and accuracy.
slide-36
SLIDE 36

Complement System

  • Complement system is a critical part of the immune system

Ricklin et al. 2007

slide-37
SLIDE 37

Classical pathway Lectin pathway Amplified response Inhibition

slide-38
SLIDE 38

Goals

  • Quantitatively understand the regulatory mechanisms of

complement system

  • How is the excessive response of the complement avoided?
  • The model:
  • Classical pathway + the lectin pathway
  • Inhibitory mechanism

 C4BP

slide-39
SLIDE 39
  • ODE Model
  • 42 Species
  • 45 Reactions

 Mass law  Michaelis-Menten kinetics

  • 92 Parameters (71 unknown)
  • DBN Construction
  • Settings

 6 intervals  100s time-step, 12600s  2.4 x 106 samples

  • Runtime

 12 hours on a cluster of 20 PCs

  • Model Calibration:
  • Training data: 4 proteins, 7 time points, 4 experimental conditions
  • Test data: Zhang et al, PLoS Pathogens, 2009

Complement System

slide-40
SLIDE 40

40

Model Calibration (parameter estimation)

slide-41
SLIDE 41

41

Model validation

 Validated the model using previous published data (Zhang et al

2009)

slide-42
SLIDE 42

42

Enhancement mechanism

 The antimicrobial response is sensitive to the pH and calcium

level

slide-43
SLIDE 43

Analysis.

 (Local and global) sensitivity analysis.  in silico experiments.

slide-44
SLIDE 44

Model predictions: The regulatory effect

  • f C4BP
  • C4BP maintains classical complement activation but

delays the maximal response time

  • But attenuates the lectin pathway activation

Classical pathway Lectin pathway

slide-45
SLIDE 45

The regulatory mechanism of C4BP

  • The major inhibitory role of C4BP is to facilitate the decay of

C3 convertase

A B D C A B C D

slide-46
SLIDE 46

Results

  • Both predictions concerning C4BP were experimentally

verified.

[PLoS Comp.Biol (2011)] [BioModels database (303.Liu)]

slide-47
SLIDE 47

Some extensions

  • Parametrized version of FF
  • Reduce errors by investing more computational time

[CMSB’11, TCBB 2012]

  • GPU implementation:
  • Significant increase in performance and scalability
  • Thrombin-dependent MLC p-pathway
  • 105 ODEs; 197 rate constants ; 164 “unknown” rate

constants.

  • (FF based approximate) probabilistic

verification method [Bioinformatics 2012]

slide-48
SLIDE 48

Current Collaborations

Immune system signaling during Multiple infections DNA damage/response pathways Chromosome co-localizations and co-regulations

Ding Jeak Ling Marie-Veronique Clement G V Shivashankar

slide-49
SLIDE 49

Conclusion

  • The DBN approximation method is useful and

efficient.

  • When does it (not) work?
  • How to relate ODEs based dynamical properties

to the DBN based ones?

  • How to extend the approximation method to

multi-mode signaling pathways?

slide-50
SLIDE 50

Acknowledgements:

David Hsu Suchee Palaniappan Liu Bing Blaise Genest Benjamin Gyori Gireedhar Venkatachalam Wang Junjie