Computational Biology, Part 20 Stochastic Modeling / Neuronal Modeling

1. Computational Biology, Part 20Stochastic Modeling / Neuronal Modeling Arvind Rao, Robert F. Murphy, Shann-Ching Chen, Justin Newberg Copyright  2004-2009. All rights reserved.

2. Stochastic Modeling in Biology

3. Stochastic Modeling in Biology • Why? A: Better resolution in species amounts • When? A: Biochemical kinetics, gene expression stochasticity in cells • How? A: SSA • Case studies: • Chemical master equation • Biochemical kinetics • Gene networks

4. Why? • Recall Chemical kinetics examples • In differential/difference equations, we examined regimes where number of molecules of reactants were always large enough to have followed mass action kinetics. • Think “Law of Large Numbers”

5. Why Are Stochastic Models Needed? • Much of the mathematical modeling of biochemical/gene networks represents gene expression deterministically • Deterministic models describe macroscopic behavior; but many cellular constituents are present in small numbers • Considerable experimental evidence indicates that significant stochastic fluctuations are present • There are many examples when deterministic models are not adequate

6. Stochastic Chemical Kinetics

9. The Chemical Master Equation

10. Challenges in the solution of CME

11. Exploiting Underlying Biology

12. Monte Carlo Simulations: Stochastic Simulation Algorithm

13. Gillespie Algorithm • Gillespie algorithm allows a discrete and stochastic simulation of a system with few reactants because every reaction is explicitly simulated. • a Gillespie realization represents a random walk that exactly represents the distribution of the Master equation. • The physical basis of the algorithm is the collision of molecules within a reaction vessel (well mixed). • all reactions within the Gillespie framework must involve at most two molecules. Reactions involving three molecules are assumed to be extremely rare and are modeled as a sequence of binary reactions.

14. Algorithm Summary • Initialization: Initialize the number of molecules in the system, reactions constants, and random number generators. • Monte Carlo Step: Generate random numbers to determine 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. • Update: Increase the time step by the randomly generated time in Step 1. Update the molecule count based on the reaction that occurred. • Iterate: Go back to Step 1 unless the number of reactants is zero or the simulation time has been exceeded.

16. SSA • Advantages • Low memory requirement • Computation is not O(exp(N)) • Disadvantages • Convergence is slow • Little insight

17. References • Daniel T. Gillespie (1977). "Exact Stochastic Simulation of Coupled Chemical Reactions". The Journal of Physical Chemistry81 (25): 2340-2361. doi:10.1021/j100540a008. • Daniel T. Gillespie (2007). “Stochastic Simulation of Chemical Kinetics". Annu. Rev. Phys. Chem. 2007.58:35-55. 10.1146/annurev.physchem.58.032806.104637 • D.Wilkinson (2009), “Stochastic modelling for quantitative description of heterogeneous biological systems.” Nature Reviews Genet. Feb 2009; 10(2):122-33. • Slides adapted from: http://www.cds.caltech.edu/~murray/wiki/images/d/d9/Khammash_master-15aug06.pdf • Gillespie: http://en.wikipedia.org/wiki/Gillespie_algorithm

18. Neuronal Modeling: The Hodgkin Huxley Equations

19. Basic neurophysiology • An imbalance of charge across a membrane is called a membrane potential • The major contribution to membrane potential in animal cells comes from imbalances in small ions (e.g., Na, K) • The maintenance of this imbalance is an active process carried out by ion pumps

20. Basic neurophysiology • Ion pumps require energy (ATP) to carry ions across a membrane up a concentration gradient (they generate a potential) • Ion channels allow ions to flow across a membrane down a concentration gradient (they dissipate a potential)

21. Basic neurophysiology • Example electrochemical gradients (left) • Example ion channel (right) Johnston & Wu, Foundations of Cellular Neurophysiology, 5th ed.

22. Basic neurophysiology • The cytoplasm of most cells (including neurons) has an excess of negative ions over positive ions (due to active pumping of sodium ions out of the cell) • By convention this is referred to as a negative membrane potential (inside minus outside) • Typical resting potential is -50 mV

23. Basic neuro-physiology • An idealized neuron consists of • soma or cell body • contains nucleus and performs metabolic functions • dendrites • receive signals from other neurons through synapses • axon • propagates signal away from soma • terminal branches • form synapses with other neurons

24. Basic neurophysiology • The junction between the soma and the axon is called the axonhillock • The soma sums (“integrates”) currents (“inputs”) from the dendrites • When the received currents result in a sufficient change in the membrane potential, a rapid depolarization is initiated in the axon hillock

25. Basic neuro-physiology • Electrical signals regulate local calcium concentrations • Synaptic vesicles fuse with the axon membrane, and neurotransmitters are released into the space between axon and dendrite in a process mediated by calcium ions • Binding of neurotransmitters to dendrite causes influx of sodium ions that diffuse into soma

26. Basic electrophysiology • A cell is said to be electrically polarized when it has a non-zero membrane potential • A dissipation (partial or total) of the membrane potential is referred to as a depolarization, while restoration of the resting potential is termed repolarization

27. Action potential in neurons • 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 Voltage difference (inside – outside) Potassium (K+) inside time

28. Action potential is linked to ion channel conductances G is channel conductance. High conductance allows for ions to pass through channel easier.

29. The Hodgkin-Huxley model • Based on electrophysiological measurements of giant squid axon • Empirical model that predicts experimental data with very high degree of accuracy • Provides insight into mechanism of action potential http://www.mun.ca/biology/desmid/brian/BIOL2060/BIOL2060-13/1310.jpg

30. The Hodgkin-Huxley model • Define • v(t)  voltage across the membrane at time t • q(t) net charge inside the neuron at t • I(t) current of positive ions into neuron at t • g(v) conductance of membrane at voltage v • Ccapacitance of the membrane • Subscripts Na, K and L used to denote specific currents or conductances (L=“other”) (INa , IK , IL ) (gNa , gK , gL )

31. The Hodgkin-Huxley model Note: Conductance is 1/R, where R is resistance E indicates membrane potential, Exare equilibrium potentials Experiments show only gNa and gK vary with time when stimulus is applied

32. The Hodgkin-Huxley model • Start with equation for capacitor

33. The Hodgkin-Huxley model • Consider each ion separately and sum currents to get rate of change in charge and hence voltage

34. The Hodgkin-Huxley model • Central concept of model: Define three state variables that represent (or “control”) the opening and closing of ion channels • m controls Na channel opening • h controls Na channel closing • n controls K channel opening

35. The Hodgkin-Huxley model • Define relationship of state variables to conductances of Na and K Q: How were the powers determined? A: Smart guessing m h n

36. The Hodgkin-Huxley model • Define empirical differential equations to model behavior of each gate

37. The Hodgkin-Huxley model • Define empirical differential equations to model behavior of each gate

38. The Hodgkin-Huxley model • Define empirical differential equations to model behavior of each gate

39. The Hodgkin-Huxley model • Gives set of four coupled, non-linear, ordinary differential equations • Must be integrated numerically • Constants (g in mmho/cm2 and v in mV)

40. Hodgkin-Huxley gates

41. Interactive demonstration • (Integration of Hodgkin-Huxley equations using Maple)

42. Interactive demonstration > Ena:=55: Ek:=-82: El:= -59: gkbar:=24.34: gnabar:=70.7: > gl:=0.3: vrest:=-69: cm:=0.001: > 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)+ 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):

43. 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(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}, {V(t),n(t),m(t),h(t)},type=numeric, 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(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}, {V(t),n(t),m(t),h(t)},type=numeric); > with(plots):

44. Interactive demonstration > 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: > 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(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}, {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});

45. Interactive demonstration > L:=odeplot(sol20,[V(t),n(t)],0..0.02,numpoints=400, 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}, {x(t),y(t)},type=numeric, output=listprocedure); > 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});

46. Virtual Cell - Hodgkin-Huxley • Versions of the models in “Computational cell biology” by Fall et al have been implemented in Virtual Cell • These are available as Public models • The Hodgkin-Huxley Model is a scientific model that describes how action potentials in neurons are initiated and propagated • Within Virtual Cell, use File/Open/Biomodel • Then open Shared Models/CompCell/Hodgkin-Huxley

47. 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"