LecturesPart20 - Computational Biology, Part 20 Stochastic...

Info iconThis preview shows page 1. Sign up to view the full content.

View Full Document Right Arrow Icon
This is the end of the preview. Sign up to access the rest of the document.

Unformatted text preview: Computational Biology, Part 20 Stochastic Modeling / Neuronal Modeling Arvind Rao, Robert F. Murphy, ShannChing Chen, Justin Newberg Copyright © 2004-2009. 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 2340-2361. doi:10.1021/j100540a008. doi Daniel T. Gillespie (2007). “Stochastic Simulation of Chemical Daniel Kinetics". Annu. Rev. Phys. Chem. 2007.58:35-55. 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):122-33. 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 non-zero 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.mcgraw-hill.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 Hodgkin-Huxley 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 L2060-13/1310.jpg The Hodgkin-Huxley 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 Hodgkin-Huxley 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 Hodgkin-Huxley model s Start with equation for capacitor v(t ) = q(t ) C The Hodgkin-Huxley 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 Hodgkin-Huxley 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 Hodgkin-Huxley 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 Hodgkin-Huxley 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 Hodgkin-Huxley 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 Hodgkin-Huxley 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 Hodgkin-Huxley model s Gives set of four coupled, non-linear, 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 Hodgkin-Huxley 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 Hodgkin-Huxley 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-(v-vrest))/(exp(0.1*(10-(v-vrest)))1): > betan:=v-> 0.125*exp(-(v-vrest)/80): > alpham:=v-> 0.1*(25-(v-vrest))/(exp(0.1*(25-(v-vrest)))-1): > betam:=v-> 4*exp(-(v-vrest)/18): > alphah:=v->0.07*exp(-0.05*(v-vrest)): > betah:=v->1/(exp(0.1*(30-(v-vrest)))+1): > pulse:=t->-20*(Heaviside(t-.001)-Heaviside(t-.002)): > rhsV:=(t,V,n,m,h)->-(gnabar*m^3*h*(V-Ena) + > gkbar*n^4*(V-Ek) + gl*(V-El) gkbar*n^4*(V-Ek) + pulse(t))/cm: > rhsn:=(t,V,n,m,h)-> 1000*(alphan(V)*(1-n) - betan(V)*n): > rhsm:=(t,V,n,m,h)-> 1000*(alpham(V)*(1-m) - betam(V)*m): > rhsh:=(t,V,n,m,h)-> 1000*(alphah(V)*(1-h) - 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*(V-Ena) + gkbar*n^4*(V-Ek) + gl*(V-El)+ pulse(t))/cm: gkbar*n^4*(V-Ek) > 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)->x-x^3/3-y; > rhsy:=(t,x,y)->c*(x+a-b*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,V-V^3/3,V=-2.5..2.2]}): > plots[display]({J,K}); Virtual Cell - Hodgkin-Huxley 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 Hodgkin-Huxley Model is a scientific The Hodgkin-Huxley 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/Hodgkin-Huxley Models/CompCell/Hodgkin-Huxley Compartments • Extracellular • Plasma Membrane • Cytosol Biochemical Species • Potassium • Sodium • Potassium Channel Inactivation Gate-closed "n_c" • Potassium Channel Inactivation Gate-open "n_o" • Sodium Channel Inactivation Gate-closed "h_c" • Sodium Channel Inactivation Gate-open "h_o" • Sodium Channel Activation Gate-open "m_o" • Sodium Channel Activation Gate-closed "m_c" m_o: Sodium Channel Activation Gate-open ...
View Full Document

Ask a homework question - tutors are online