Transcript Document

Gillespie’s Framework
Track precise (integer) quantities of molecular species.
“States”
A
Reactions
B
C
S1 4
7
5
R1
S2 2
6
8
R2
S3 22
0
997
R3




k1
2A  B  3 C
k2
B  2C 
A  C  2B
A reaction transforms one state into another:
S2
e.g., S1 
R
1
k3
3A
Stochastic
Simulation
S1 = [5, 5, 5] 0
R1 R2
R3
Ri
ki
ni ,1 X i ,1  ni , 2 X i , 2   
Choose the next reaction
according to:
i
Pr( Ri ) 
 j
j
where
 X 1  X 2 

 i  k i  
 n1  n2 
Stochastic
Simulation
S1 = [5, 5, 5] 0
R1 R2
R3
Ri
ki
ni ,1 X i ,1  ni , 2 X i , 2   
Choose the time of the next
reaction according to:
Pr(t  t0 )  
t0
 0
 je
j





j

 j 


d
S1 = [5, 5, 5] 0
R1 R2
Stochastic
Simulation
R3
See D. Gillespie,
“Exact Stochastic Simulation of
Coupled Chemical Reactions”,
J. Phys. Chem. 1977
S1 = [5, 5, 5] 0
Stochastic
Simulation
Choose R3 and t = 3 seconds.
R1 R2
R3
S2 = [4, 7, 4] 3
Choose R1 and t = 1 seconds.
S3 = [2, 6, 7] 4
Choose R3 and t = 2 seconds.
S4 = [1, 8, 6] 6
Choose R2 and t = 1 seconds.
S1 = [5, 5, 5] 0
Stochastic
Simulation
Choose R3 and t = 3 seconds.
S2 = [4, 7, 4] 37
Choose R1 and t = 1 seconds.
S3 = [2, 6, 7] 4
Choose R3 and t = 2 seconds.
S4 = [1, 8, 6] 6
Choose R2 and t = 1 seconds.
S1 = [5, 5, 5] 0
Stochastic
Simulation
S2 = [4, 7, 4] 7
S3 = [2, 6, 7] 4
S4 = [1, 8, 6] 6
Stochastic simulation
spends most of its time
(e.g. 99%) looping
Is Looping Necessary?
A+N
Reactions
B+N
C
1
2
3
B
C
A + 2N
0.001
B
X
0.002
C
Y
If we begin with 1 molecule of A and 2 molecules of N:
• what is the probability that we get a molecule of X?
• what is the probability that we get a molecule of Y?
(assuming that we wait as long as it takes)
Is Looping Necessary?
A+N
Reactions
B+N
C
Trial 2
1
SA
Cycle 1000 times ...
Now repeat 500 times ...
1
2
3
B
0.001
B
C
X
0.002
C
A + 2N
SB
Y
SC
SX
SY
Is Looping Necessary?
1
A+N
Reactions
2
B+N
3
C
SA
1
Try computing probabilities!
B
0.001
B
C
X
0.002
C
A + 2N
SB
p
Y
SC
q
1 p
1 q
SX
SY
1 p
p (1  q )
Is Looping Necessary?
Pr( S X ) 
Conclude:
1 p
1 p

(1  p )  p (1  q ) 1  pq
p (1  q )
p (1  q )
Pr( SY ) 

(1  p )  p (1  q ) 1  pq
SA
1
Try computing probabilities!
SB
p
SC
q
1 p
1 q
SX
SY
1 p
p (1  q )
Is Looping Necessary?
Pr( S X ) 
Conclude:
1 p
1 p

(1  p )  p (1  q ) 1  pq
p (1  q )
p (1  q )
Pr( SY ) 

(1  p )  p (1  q ) 1  pq
SA
Break the cycle
SB
SC
SX
SY
Breaking Cycles
Whenever a cycle is encountered during stochastic simulation:
•
•
•
Compute the probability from the entry point to each exit point.
Create new transitions to the exit points with these probabilities.
Break the cycle at the entry point.
Breaking Cycles
Whenever a cycle is encountered during stochastic simulation:
•
•
•
Compute the probability from the entry point to each exit point.
Create new transitions to the exit points with these probabilities.
Break the cycle at the entry point.
Breaking Cycles
Caveats:
• We must keep a history of states visited.
(However, we should be caching state information anyway...)
• We’ve lost the concept of time.
(This can be rectified with a bit of math....)
Randomness
Pseudo-random
numbers needed:
probabilities
R1
R2
R3
R4
1
16
1
1
16
3
8
4
0
generate a random number: 0.07123  choose R2
1
Randomness
Pseudo-random
numbers needed:
R1
R2
R3
R4
probabilities
0
chooseR4R2
generate a random number: 0.07123
0.8973 choose
1
Randomness
Pseudo-random
numbers needed:
R1
R2
R3
R4
probabilities
• Generating random numbers is time consuming.
• If variance in probabilities is large, accuracy is wasted.
0
1
generate a random number: 0.8973  choose R4
Event Leaping
Explore high probability
events further.
1
16
1
8
1
16
3
4
Along each path, probabilities
are multiplicative.
3
16
3
16
3
16
3
16
Event Leaping
Explore high probability
events further.
1
16
1
8
1
16
Along each path, probabilities
are multiplicative.
1
When paths merge,
32
probabilities are additive.
7
3
16
32
3
16
7
3
16
32
7
3
16
32
Event Leaping
Exploreon
Based
high
a single
probability
randomfurther.
events
number, leap
directly to the boundary
of explored region.
1
16
1
16
1
When paths merge,
32
probabilities are additive.
3
16
7
32
7
32
7
32
Event Leaping
• We need far fewer random numbers (e.g., factor of
10 reduction).
• We cache probability calculations.
• If we return to any portion of the region already
visited, we immediately jump to the boundary of it.
State Space Analysis
start
[3, 3, 3]
Characterize Evolution
e.g., identify articulation points
p1
p2
[2, 2, 6]
p4
[1, 1, 9]
p10
[0, 0, 12]
p1 p 4 p10
p3
[2, 6, 1]
[5, 1, 2]
p6
p7
p5
[1, 5, 4]
p8
[4, 4, 0]
p9
[4, 0, 5]
p12
p11
p13
[3, 3, 3]
p11 ( p 1 p5  p 2 p 6 )  p12 ( p 2 p 7  p3 p8 )  p3 p9 p13
State Space Analysis
Characterize Evolution
e.g., identify “articulation” points
[2, 2, 6]
Pr(C1) > 0.99
start
[3, 3, 3]
[2, 6, 1]
[5, 1, 2]
Pr(C2) > 0.99