This preview shows page 1. Sign up to view the full content.
Unformatted text preview: Computational Biology, Part 20
Stochastic Modeling / Neuronal
Modeling
Arvind Rao, Robert F. Murphy, ShannChing Chen, Justin Newberg
Copyright © 20042009.
Copyright
All rights reserved. Stochastic Modeling in Biology Stochastic Modeling in Biology
s Why? A: Better resolution in species
Why?
resolution
amounts
amounts
s When? A: Biochemical kinetics, gene
When?
kinetics gene
expression stochasticity in cells
expression
s How? A: SSA
How?
SSA
s Case studies:
x Chemical master equation
3 Biochemical kinetics
3 Gene networks s Why?
s Recall Chemical kinetics examples
1. In differential/difference equations, we
In
examined regimes where number of
molecules of reactants were always large
enough to have followed mass action
kinetics.
kinetics. 2. Think “Law of Large Numbers” Why Are Stochastic Models
Needed?
• Much of the mathematical modeling of
biochemical/gene networks represents gene expression
deterministically
deterministically
• Deterministic models describe macroscopic behavior;
but many cellular constituents are present in small
numbers
numbers
• Considerable experimental evidence indicates that
significant stochastic fluctuations are present
significant
• There are many examples when deterministic models
are not adequate
not Stochastic Chemical Kinetics The Chemical Master Equation Challenges in the solution of
CME Exploiting Underlying Biology Monte Carlo Simulations:
Stochastic Simulation Algorithm Gillespie Algorithm
s s s s Gillespie algorithm allows a discrete and stochastic
Gillespie
simulation of a system with few reactants because every
reaction is explicitly simulated.
a Gillespie realization represents a random walk that exactly
Gillespie
represents the distribution of the Master equation.
represents
The physical basis of the algorithm is the collision of
The
molecules within a reaction vessel (well mixed).
all reactions within the Gillespie framework must involve at
all
most two molecules. Reactions involving three molecules are
assumed to be extremely rare and are modeled as a sequence
of binary reactions. Algorithm Summary
1.
1. 2. 3. 4. Initialization: Initialize the number of molecules in the
Initialize
system, reactions constants, and random number generators.
system,
Monte Carlo Step: Generate random numbers to determine
Generate
the next reaction to occur as well as the time interval. The
probability of a given reaction to be chosen is proportional
to the number of substrate molecules.
to
Update: Increase the time step by the randomly generated
Increase
time in Step 1. Update the molecule count based on the
reaction that occurred.
reaction
Iterate: Go back to Step 1 unless the number of reactants is
Go
zero or the simulation time has been exceeded.
zero SSA
s Advantages
x Low memory requirement
x Computation is not O(exp(N)) s Disadvantages
x Convergence is slow
x Little insight References
s s s s s Daniel T. Gillespie (1977). "Exact Stochastic Simulation of Coupled
Daniel
Chemical Reactions". The Journal of Physical Chemistry 81 (25):
The
81
23402361. doi:10.1021/j100540a008.
doi
Daniel T. Gillespie (2007). “Stochastic Simulation of Chemical
Daniel
Kinetics". Annu. Rev. Phys. Chem. 2007.58:3555.
10.1146/annurev.physchem.58.032806.104637
10.1146/annurev.physchem.58.032806.104637
D.Wilkinson (2009), “Stochastic modelling for quantitative description
D.Wilkinson
of heterogeneous biological systems.” Nature Reviews Genet. Feb
2009; 10(2):12233.
Slides adapted from:
Slides
http://www.cds.caltech.edu/~murray/wiki/images/d/d9/Khammash_masterhttp://www.cds.caltech.edu/~murray/wiki/images/d/d9/Khammash_masterGillespie: http://en.wikipedia.org/wiki/Gillespie_algorithm
Gillespie: http://en.wikipedia.org/wiki/ Neuronal Modeling: The Hodgkin
Huxley Equations Basic neurophysiology
s An imbalance of charge across a membrane
An
is called a membrane potential
membrane
s The major contribution to membrane
The
potential in animal cells comes from
imbalances in small ions (e.g., Na, K)
imbalances
s The maintenance of this imbalance is an
The
active process carried out by ion pumps
active Basic neurophysiology
s Ion pumps require energy (ATP) to carry
Ion
ions across a membrane up a concentration
up
gradient (they generate a potential)
potential)
s Ion channels allow ions to flow across a
Ion
membrane down a concentration gradient
down
(they dissipate a potential)
dissipate Basic neurophysiology
s Example electrochemical gradients (left)
s Example ion channel (right)
Johnston & Wu, Foundations of Cellular
Johnston
Neurophysiology, 5th ed.
Neurophysiology Basic neurophysiology
s The cytoplasm of most cells (including
The
neurons) has an excess of negative ions over
positive ions (due to active pumping of
sodium ions out of the cell)
sodium
s By convention this is referred to as a
By
negative membrane potential (inside
minus outside)
minus
s Typical resting potential is 50 mV
Typical
is Basic neurophysiology
s An idealized neuron consists of
An
neuron
x soma or cell body
cell
3 contains nucleus and performs metabolic functions x dendrites
3 receive signals from other neurons through synapses
receive
synapses x axon
3 propagates signal away from soma x terminal branches
3 form synapses with other neurons
form synapses Basic neurophysiology
s The junction between the soma and the axon
The
is called the axon hillock
axon hillock
s The soma sums (“integrates”) currents
The
(“inputs”) from the dendrites
(“inputs”)
s When the received currents result in a
When
sufficient change in the membrane potential,
a rapid depolarization is initiated in the axon
hillock
hillock Basic neurophysiology
s s s Electrical signals regulate local calcium
Electrical
calcium
concentrations
concentrations
Synaptic vesicles fuse with the axon membrane,
and neurotransmitters are released into the space
neurotransmitters
between axon and dendrite in a process mediated
by calcium ions
by
Binding of neurotransmitters to dendrite causes
Binding
influx of sodium ions that diffuse into soma Basic electrophysiology
s A cell is said to be electrically polarized
cell
polarized
when it has a nonzero membrane potential
when
s A dissipation (partial or total) of the
dissipation
membrane potential is referred to as a
depolarization, while restoration of the
depolarization while
resting potential is termed repolarization
repolarization Action potential in neurons
s Voltage difference (inside – outside) s http://highered.mcgrawhill.com/olc/dl/120107/bio_d.swf
http://bcs.whfreeman.com/thelifewire/content/chp44/4402002.html outside
Sodium (Na+)
Cell Membrane inside time Potassium (K+) Action potential is linked to ion
channel conductances
Sti mu l u s (u A ) 150
100 G is channel
is
conductance.
conductance. 50
0
60 Vo l tag e (mV) 40 High conductance
High
allows for ions to
pass through
channel easier.
channel 20
0
20
40
60 C o n d u ctan ce (mS/ cm2) 80
40
G(Na) 30 G(K) 20
10
0
0 2 4 Time (ms) 6 8 10 The HodgkinHuxley model
s s s Based on
Based
electrophysiological
measurements of giant
squid axon
squid
Empirical model that
Empirical
predicts experimental data
with very high degree of
accuracy
accuracy
Provides insight into
Provides
mechanism of action
potential
potential
http://www.mun.ca/biology/desmid/brian/BIOL2060/BIO
L206013/1310.jpg The HodgkinHuxley model
s Define
x v(t) ≡
v(t) voltage across the membrane at time t
x q(t) ≡ net charge inside the neuron at t
q(t)
x I(t) ≡ current of positive ions into neuron at t
I(t)
x g(v) ≡ conductance of membrane at voltage v
g(v)
x C ≡ capacitance of the membrane
capacitance
x Subscripts Na, K and L used to denote specific
Subscripts
currents or conductances (L=“other”)
currents
(INa , IK , IL ) (gNa , gK , gL ) The HodgkinHuxley model
Note:
Conductance is
Conductance
1/R, where R is
1/R where
resistance
resistance
E indicates
membrane
potential, Ex are
equilibrium
potentials
potentials
Experiments
Experiments
show only gNa
and gK vary with
time when
stimulus is
applied
applied The HodgkinHuxley model
s Start with equation for capacitor
v(t ) = q(t )
C The HodgkinHuxley model
s Consider each ion separately and sum
Consider
currents to get rate of change in charge and
hence voltage
hence
dq
= −( I Na + I K + I L )
dt
I Na = g Na (v − v Na )
I K = g K (v − v K )
I L = g L (v − v L )
dv −1
[ g Na (v)(v − vNa ) + g K (v)(v − vK ) + g L (v − vL )]
=
dt
C The HodgkinHuxley model
s Central concept of model: Define three state
Central
variables that represent (or “control”) the
opening and closing of ion channels
opening
x m controls Na channel opening
x h controls Na channel closing
x n controls K channel opening The HodgkinHuxley model
s Define relationship of state variables to
Define
conductances of Na and K
conductances
Q: How were
Q:
the powers
g Na = g Na m 3 h
gK = gK n determined? 4 A: Smart
guessing
guessing 0 ≤ m, n, h ≤ 1 m h n The HodgkinHuxley model
s Define empirical differential equations to
Define
model behavior of each gate
model dn
= α n (v )(1 − n) − βn (v )n
dt
0.01(v + 10)
α n (v) =
( v +10 ) /10
(e
− 1) βn (v) = 0.125e v / 80 The HodgkinHuxley model
s Define empirical differential equations to
Define
model behavior of each gate
model dm
= α m (v )(1 − m) − βm (v )m
dt
0.1(v + 25)
α m (v) =
( v + 25 ) /10
(e
− 1) βm (v) = 4e v /18 The HodgkinHuxley model
s Define empirical differential equations to
Define
model behavior of each gate
model dh
= α h (v )(1 − h ) − β h (v ) h
dt
v / 20
α h (v ) = 0.07e
β h (v ) = 1 ( v + 30) /10
(e
+ 1) The HodgkinHuxley model
s Gives set of four coupled, nonlinear,
Gives
ordinary differential equations
ordinary
s Must be integrated numerically
s Constants (g in mmho/cm2 and v in mV)
in g Na = 120 vNa = −115
g K = 36
g L = 0 .3 vK = 12
vL = −10.5989 HodgkinHuxley gates
Sti mu l u s (u A ) 150
100
50
0
60 150 20
0 Sti mu l u s (u A ) Vo l tag e (mV) 40 100 20
40
60
80 G(Na) 30 G(K) 20 40 0
2 4 Time (ms) 6 8 10 Vo l tag e (mV) 0 0
60 10 20
0
20
40
60
80
1.0 Gate p ar am val u e C o n d u ctan ce (mS/ cm2) 40 50 m gate (Na) 0.8 h gate (Na) 0.6 n gate (K) 0.4
0.2
0.0
0 2 4 Time (ms) 6 8 10 Interactive demonstration
s (Integration of HodgkinHuxley equations
(Integration
using Maple)
using Interactive demonstration
> Ena:=55: Ek:=82: El:= 59: gkbar:=24.34: gnabar:=70.7:
Ena:=55:
> gl:=0.3: vrest:=69: cm:=0.001:
gl:=0.3:
> alphan:=v> 0.01*(10(vvrest))/(exp(0.1*(10(vvrest)))1):
> betan:=v> 0.125*exp((vvrest)/80):
> alpham:=v> 0.1*(25(vvrest))/(exp(0.1*(25(vvrest)))1):
> betam:=v> 4*exp((vvrest)/18):
> alphah:=v>0.07*exp(0.05*(vvrest)):
> betah:=v>1/(exp(0.1*(30(vvrest)))+1):
> pulse:=t>20*(Heaviside(t.001)Heaviside(t.002)):
> rhsV:=(t,V,n,m,h)>(gnabar*m^3*h*(VEna) +
>
gkbar*n^4*(VEk) + gl*(VEl)
gkbar*n^4*(VEk)
+ pulse(t))/cm:
> rhsn:=(t,V,n,m,h)> 1000*(alphan(V)*(1n)  betan(V)*n):
> rhsm:=(t,V,n,m,h)> 1000*(alpham(V)*(1m)  betam(V)*m):
> rhsh:=(t,V,n,m,h)> 1000*(alphah(V)*(1h)  betah(V)*h): Interactive demonstration
> inits:=V(0)=vrest,n(0)=0.315,m(0)=0.042, h(0)=0.608;
> sol:=dsolve({diff(V(t),t)=rhsV(t,V(t),n(t),m(t),h(t)),
diff(n(t),t)=rhsn(t,V(t),n(t),m(t),h(t)),
diff(n(t),t)=rhsn(t,V(t),n(t),m(t),h(t)),
diff(m(t),t)=rhsm(t,V(t),n(t),m(t),h(t)),
diff(m(t),t)=rhsm(t,V(t),n(t),m(t),h(t)),
diff(h(t),t)=rhsh(t,V(t),n(t),m(t),h(t)),inits},
diff(h(t),t)=rhsh(t,V(t),n(t),m(t),h(t)),inits},
{V(t),n(t),m(t),h(t)},type=numeric,
output=listprocedure);
output=listprocedure);
> Vs:=subs(sol,V(t));
> plot(Vs,0..0.02);
> sol20:=dsolve({diff(V(t),t)=rhsV(t,V(t),n(t),m(t),h(t)),
diff(n(t),t)=rhsn(t,V(t),n(t),m(t),h(t)),
diff(n(t),t)=rhsn(t,V(t),n(t),m(t),h(t)),
diff(m(t),t)=rhsm(t,V(t),n(t),m(t),h(t)),
diff(m(t),t)=rhsm(t,V(t),n(t),m(t),h(t)),
diff(h(t),t)=rhsh(t,V(t),n(t),m(t),h(t)),inits},
diff(h(t),t)=rhsh(t,V(t),n(t),m(t),h(t)),inits},
{V(t),n(t),m(t),h(t)},type=numeric);
{V(t),n(t),m(t),h(t)},type=numeric);
> with(plots): Interactive demonstration
> J:=odeplot(sol20,[V(t),n(t)],0..0.02):
J:=odeplot(sol20,[V(t),n(t)],0..0.02): > display({J});
> pulse:=t>2*(Heaviside(t.001)Heaviside(t.002)):
> rhsV:=(t,V,n,m,h)>(gnabar*m^3*h*(VEna) +
gkbar*n^4*(VEk) + gl*(VEl)+ pulse(t))/cm:
gkbar*n^4*(VEk)
> sol2:=dsolve({diff(V(t),t)=rhsV(t,V(t),n(t),m(t),h(t)),
diff(n(t),t)=rhsn(t,V(t),n(t),m(t),h(t)),
diff(n(t),t)=rhsn(t,V(t),n(t),m(t),h(t)),
diff(m(t),t)=rhsm(t,V(t),n(t),m(t),h(t)),
diff(m(t),t)=rhsm(t,V(t),n(t),m(t),h(t)),
diff(h(t),t)=rhsh(t,V(t),n(t),m(t),h(t)),inits},
diff(h(t),t)=rhsh(t,V(t),n(t),m(t),h(t)),inits},
{V(t),n(t),m(t),h(t)},type=numeric);
{V(t),n(t),m(t),h(t)},type=numeric);
> K:=odeplot(sol2,[V(t),n(t)],0..0.02,color=green):
> display({J,K}); Interactive demonstration
> L:=odeplot(sol20,[V(t),n(t)],0..0.02,numpoints=400,
color=blue):
color=blue):
> display({J,L});
> odeplot(sol20,[V(t),m(t)],0..0.02,numpoints=400);
> odeplot(sol20,[V(t),h(t)],0..0.02,numpoints=400);
> odeplot(sol20,[m(t),h(t)],0..0.02,numpoints=400);
> a:=0.7; b:=0.8; c:=0.08;
> rhsx:=(t,x,y)>xx^3/3y;
> rhsy:=(t,x,y)>c*(x+ab*y);
> sol2:=dsolve({diff(x(t),t)=rhsx(t,x(t),y(t)),
diff(y(t),t)=rhsy(t,x(t),y(t)),x(0)=0,y(0)=1},
diff(y(t),t)=rhsy(t,x(t),y(t)),x(0)=0,y(0)=1},
{x(t),y(t)},type=numeric, output=listprocedure);
{x(t),y(t)},type=numeric,
> xs:=subs(sol2,x(t)); ys:=subs(sol2,y(t));
> K:=plot([xs,ys,0..200],x=3..3,y=2..2,color=blue):
> J:=plot({[V,(V+a)/b,V=2.5..1.5],[V,VV^3/3,V=2.5..2.2]}):
> plots[display]({J,K}); Virtual Cell  HodgkinHuxley
s Versions of the models in “Computational
Versions
cell biology” by Fall et al have been
implemented in Virtual Cell
implemented
s These are available as Public models
s The HodgkinHuxley Model is a scientific
The HodgkinHuxley
model that describes how action potentials in
neurons are initiated and propagated
neurons
s Within Virtual Cell, use File/Open/Biomodel
s Then open Shared
Then
Models/CompCell/HodgkinHuxley
Models/CompCell/HodgkinHuxley Compartments
• Extracellular
• Plasma Membrane
• Cytosol Biochemical Species
• Potassium
• Sodium
• Potassium Channel Inactivation Gateclosed "n_c"
• Potassium Channel Inactivation Gateopen "n_o"
• Sodium Channel Inactivation Gateclosed "h_c"
• Sodium Channel Inactivation Gateopen "h_o"
• Sodium Channel Activation Gateopen "m_o"
• Sodium Channel Activation Gateclosed "m_c" m_o: Sodium Channel Activation Gateopen ...
View Full
Document
 Fall '08
 Staff
 Computational Biology

Click to edit the document details