Anali PAZU - Letnik 1, leto 2011, številka 2 122 Simulacija akcijskega potenciala v Hodgkin-Huxleyevem modelu Simulation oftheactionpotential in the Hodgkin-Huxley model Rudolf Pušenjak, Maks Oblak Univerza v Mariboru, Fakulteta za logistiko / Mariborska cesta 17, SI-3000 Celje E-mail: rudi.pusenjak@teleing.com Povzetek: Vzdražljive (ekscitabilne) celice se imenujejo celice, ki so sposobne proizvajati akcijske potenciale. Akcijski potenciali nastajajojo v nevronih in mišicah v obliki signalov. Signalni mehanizem povzroča možgansko aktivnost oziroma krčenje mišic. Lipidni plasti celične membrane se obnašata kot kondenzator, spremenljive prevodnosti membrane za posamezne ione pa so posledica ionskih kanalov, ki se lahko odpirajo ali zapirajo in s tem bolje ali slabše prevajajo. Na tej osnovi lahko celično membrano modeliramo kot električno vezje. V članku je predstavljena izpeljava enačb Hodgkin-Huxleyevega modela celične membrane s pomočjo osnovnih zakonov elektrotehnike in kinetike prvega reda. Dobljeni matematični model celične membrane je uporabljen za simulacijo akcijskega potenciala s simboličnim programskim orodjem Mathematica ® . Ključne besede: Hodgkin-Huxleyev model ; celična membrane ; simulacija akcijskega potencial Summary: Cells are called excitable, when produce action potentials. Action potentials arise in neurons and muscle cells in the form of signals. Signaling mechanism causes the brain activity or muscular contraction. Lipid layers of membrane cell act as a capacitor and variable membrane conductances are consequence of ion channels, which can open or close. On this basis, the cell membrane can be modelled as an electrical circuit. In the paper, equations of the Hodgkin-Huxley cell membrane model are derived by using fundamental electrical laws and first order kinetics. The gained mathematical model of the cell membrane is used in the simulation of the action potential by symbolic programming tool such as Mathematica ® . Key words: Hodgkin-Huxley model; cell membrane; simulation of action potential 1. Uvod Fiziologija nastanka in širjenja akcijskih potencialov v nevronskih in mišičnih celicah je bila predmet raziskav v vsem preteklem stoletju. Najpomembnejši znanstveni dosežek na področju teh raziskav predstavlja Hodgkin- Huxleyev kvantitativni matematični model [1], podprt z eksperimentalnimi rezultati, za katerega sta avtorja leta 1963 prejela Nobelovo nagrado v fiziologiji in medicini. Namen tega članka je predstaviti osnovne zamisli modela in njegove uporabe v simulaciji širjenja akcijskega potenciala s simboličnim jezikom programskega orodja Mathematica ® . 2. Hodgkin-Huxleyev model celične membrane Prevodnostni model je v elektrofiziologiji najpreprostejša predstava vzdražljivih celic, kakršni so nevroni oziroma celice srčne mišice. Lipidni plasti celične membrane se vedeta kot kondenzator, skozi ionske kanale membrane, ki se lahko odpirajo ali zapirajo, pa lahko prehajajo posamezni ioni (slika 1.a). Prehajanju ionov skozi kanale ustrezajo spremenljive prevodnosti, ki so odvisne od stopnje odprtosti posameznega ionskega kanala. Celotni Hodgkin- Huxleyev model celične membrane lahko predstavimo kot električno vezje (slika 1.b), v katerem teče membranski električni tok I(t) kot vsota kapacitivnega toka preko membranske kapacitivnosti C, kadar se Anali PAZU - Letnik 1, leto 2011, številka 2 123 kapacitivnost membrane spremenljivo naelektri v odvisnosti od časa in ionskih tokov, ki tečejo preko paralelno vezanih spremenljivih prevodnosti. Ionske tokove predstavljajo tokovi natrijevih (I Na ) in kalijevih ionov (I K ) ter razmeroma majhnega deleža klorovih, kalcijevih in drugih ionov, ki so združeni v tako imenovanem prepustnem (´´leakage´´) ionskem toku I L . Ionske tokove poganjajo gonilne napetosti, ki jih v vsaki posamezni paralelni veji predstavlja razlika med membranskim potencialom V in ravnotežnimi potenciali posameznih ionov E Na , E K in E L , določenimi z Nernstovo enačbo. Paralelne veje električnega vezja obravnavamo kot kanale, ki jih lahko s pomočjo spremenljivih prevodnosti g Na , g K in g L odpiramo ali zapiramo. Mehanizem odpiranja in zapiranja ionskih kanalov določajo kinetične enačbe prvega reda, s čemer postane Hodkin-Huxleyev model nelinearni dinamični sistem, ki omogoča analizo številnih pojavov nelinearnih nihanj kot so bifurkacije in celo pojav kaosa [3],[5],[6]. Slika 1. Hodgkin-Huxleyev model celične membrane Neposredno iz slike 1.b lahko s pomočjo osnovnih zakonitosti elektrotehnike napišemo naslednjo Hodgkin- Huxleyevo enačbo membranskega toka ( ) ( ) ( ) ( ) ( ) ( ) d d d d Na K L Na Na K K L L V t I t C I I I t V t C g V t E g V t E g V t E t = + + + = + − + − + −             , (1) ki potrjuje, da je model celične membrane dinamični sistem. Na osnovi eksperimentov sklepamo, da model lahko opisuje fiziološke lastnosti membrane, če sta prevodnosti g Na in g K funkciji časa in membranskega potenciala, medtem ko smemo E Na , E K , E L , C in g L vzeti konstantne. V skladu s tem najprej izrazimo ravnotežne ali Nernstove potenciale s pomočjo enačbe: ( ) ln , , , ion ion ion ec RT ion z F ic E ion Na K L   = =     , (2) kjer je R univerzalna plinska konstanta, T absolutna temperatura, z ion valenčno število iona, F Faradayeva konstanta, ion ec in ion ic pa izven celična koncentracija oziroma koncentracija ionov v notranjosti celice. Membranske prevodnosti za natrijeve in kalijeve ione so odvisne od števila kanalov in njihovega stanja. Membransko prevodnost natrijevih ionov modeliramo kot produkt maksimalne prevodnosti kanala G ion z aktivacijsko funkcijo m in deaktivacijsko funkcijo h v obliki: ( ) ( ) ( ) , , , , q r ion ion g G m V t h V t ion Na K = = , (3) kjer sta q in r eksponenta, ki funkcionalno odvisnost prevodnosti natrijevih ionov najbolje prilagodita izmerjenim podatkom. Membransko prevodnost kalijevih ionov modeliramo z eno samo aktivacijsko funkcijo n in ustreznim eksponentom q. Aktivacijska in deaktivacijska funkcija kanala v enačbi (4) sta odvisni od membranskega potenciala V in časa t in ju lahko tolmačimo kot verjetnost, da je kanal v odprtem stanju. Časovna odvisnost obeh funkcij je določena s kinetiko prvega reda in je za obe funkciji enaka. V splošnem lahko obe funkciji zaznamujemo kot verjetnost p. Verjetnost p, da je kanal v odprtem stanju, ustreza količniku števila odprtih kanalov in celotnega števila kanalov v celici. Pri tem predpostavljamo, da lahko kanal krmilimo s pomočjo vrat, ki so lahko odprta ali zaprta. Vzemimo, da α p pomeni sorazmernostni faktor odprtih, β p pa sorazmernostni faktor zaprtih kanalov. Dinamiko verjetnosti p tedaj določa razlika med stopnjo odprtih ( ) 1 p p α − in stopnjo zaprtih kanalov p p β : ( ) d 1 d p p p p p t α β = − − , (4) kjer je stopnja odprtih kanalov v stacionarnem stanju enaka stopnji zaprtih kanalov ( ) 1 p p p p α β − = s C g Na g L g K E Na E K E L izvencelični medij medij v notranjosti celice I(t) I Na (t) I K (t) I L (t) I C (t) b) Anali PAZU - Letnik 1, leto 2011, številka 2 124 pripadajočo stacionarno rešitvijo p p p p p α α β ∞ + = = . Če poznamo začetno verjetnost p 0 v začetnem času t 0 , lahko izrazimo popolno rešitev enačbe (4) v obliki: ( ) ( ) ( )( ) 0 0 exp p p p t p p p t t α β ∞ ∞   = − − − + −   (5) Enačbe Hodgkin-Huxleyevega modela celične membrane tvorimo tako, da k enačbi membranskega toka (1) dodamo kinetične enačbe verjetnosti odprtega kanala (4). Namen simulacij akcijskega potenciala je zagotovo čim boljše ujemanje rezultatov z eksperimentalnimi meritvami. Da bi to dosegli, prevodnost natrijevih ionov v enačbi (3) opišemo s pomočjo aktivacijske funkcije m z eksponentom 3 q = in deaktivacijske funkcije h z eksponentom 1 r = v obliki: 3 Na Na g G m h = , (6) prevodnost kalijevih ionov pa z eno samo aktivacijsko funkcijo n z eksponentom 4 q= : 4 K K g G n = . (7) Celotni Hodgkin-Huxleyev model celične membrane, ki omogoča numerični izračun akcijskega potenciala, lahko sedaj zapišemo v obliki ( ) L L G g = : ( ) ( ) ( ) ( ) ( ) ( )( ) ( ) ( )( ) ( ) ( )( ) ( ) 3 4 d , d d , 1 , , d d , 1 , , d d , 1 , . d Na Na K K L L m m h h n n V t C I t G m h V t E G n V t E G V t E t m V t m V t m t h V t h V t h t n V t n V t n t α β α β α β  = − − − − − −                = − −     = − −    = − −   (8) kjer so sorazmernostni faktorji odprtih oziroma zaprtih kanalov naslednje funkcije akcijskega potenciala: ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 25 10 18 25 10 20 30 10 10 10 80 10 10 , , , 4 1 1 , 0.07 , , 1 , 0.1 , , 0.125 1 V t V t m m V t V t h h V t V t V t n n V t V t V t Exp Exp V t Exp V t Exp V t V t Exp Exp α β α β α β − − − − − − − −             = =        −           = =         +                 = =         −       (9) Iz enačb (8) in (9) izhaja, da je celotni Hodgkin- Huxleyev model celične membrane strogo nelinearen. Za njegovo reševanje sta na razpolago dve poti: numerična integracija enačb ali analitično reševanje z redukcijo modela na nizko-dimenzionalni sistem. Numerična integracija enačb je učinkovit način izračuna akcijskega potenciala, ionskih tokov, aktivacijskih in deaktivacijskih funkcij v odvisnosti od časa. Numerično integracijo lahko izvajamo v programskih jezikih Mathematica, Maple in Matlab, v katerih so standardne metode integracije (Runge-Kutta, Adams, Gear, itd.) že vgrajene in uporabniku ni potrebno, da bi se ukvarjal s programiranjem teh metod. Ta možnost se izkaže še posebno pomembna pri razširitvi modela širjenja akcijskega potenciala na tkiva (npr. na srčno mišico), ki jih obravnavamo kot električno sklopljene celice. S tem postane akcijski potencial časovno in prostorsko odvisen in ga modeliramo z valovno enačbo: ( ) ( ) , , , , , , m ion V x y z t s V x y z t C I vol t σ ∂   ∇⋅ ∇ = +       ∂   , (10) ki je smiselna razširitev enačbe (1) v primeru tkiva, V(x,y,z,t) je akcijski potencial v tkivu, ∇ V(x,y,z,t) je gradient akcijskega potenciala z ozirom na krajevne spremenljivke, σ je matrika prevodnosti, ∇ ·(σ∇ V) je izvornost (divergenca) tokovnega polja, ion I je celotni ionski tok, ki teče skozi celično membrano, m s vol pa je površina celice na enoto prostornine. Razen numerične integracije je pri enačbi (10) potrebno vključiti še metodo končnih elementov (MKE) za modeliranje geometrije celične strukture (npr. oblike srčne mišice z vsemi prekati, žilnim sistemom, Purkinijevimi vlakni itd.) tako da postane problem izredno kompleksen in so zmogljiva programska orodja (Mathematica, Maple, Matlab) prava izbira za reševanje. Po drugi strani numerična integracija Anali PAZU - Letnik 1, leto 2011, številka 2 125 ne omogoča matematične analize vpliva fizioloških parametrov na dinamiko modela v teoretičnem smislu. Analizo vpliva parametrov na širjenje akcijskega potenciala si olajšamo z redukcijo števila enačb Hodgkin-Huxleyevega modela, npr. z FitzHugh- Nagumovim modelom [2], ki ga sestavljata dve enačbi. Model omogoča analitične raziskave bifurkacij, limitnih ciklov in prehoda v kaos v fazni ravnini in s pomočjo perturbacijskih metod [3], še posebej s singularno perturbacijsko metodo [4], izvajamo pa jih lahko tudi numerično z uporabo programskih paketov BIFPACK, AUTO, itd. [5]. Enačbi FitzHugh-Nagumovega modela ob upoštevanju ustreznih pretvorb zapišemo v obliki: ( ) ( )( ) ( ) ( ) d , 1 , d d , , d v f v w v v v w I t w g v w v w t λ ε μ = = − − − + = = − (11) kjer je v hitra spremenljivka (akcijski potencial), w pa počasna spremenljivka (ki ustreza aktivacijski funkciji natrijevih ionov), ε,A,λ in μ so pozitivni parametri, pri čemer velja 0<λ<1. Sistem enačb (11) ima obliko nelinearnih navadnih diferencialnih enačb prvega reda, kakršno imajo npr. enačbe Lorenzovega kaotičnega atraktorja. V [6] smo razvili večstopenjsko homotopsko perturbacijsko metodo za kaotična nihanja Lorenzovega atraktorja, ki omogoča njegovo analitično reševanje in s tem analizo bifurkacij, limitnih ciklov ter analizo stabilnosti omenjenega sistema. V nadalnjih raziskavah pričakujemo, da bomo metodo zaradi podobne zgradbe enačb lahko uspešno uporabili tudi v elektrofizioloških raziskavah na področju širjenja akcijskega potenciala. 3. Simulacija akcijskega potenciala in ionskih tokov celične membrane Simulacijo akcijskega potenciala, ionskih tokov, aktivacijskih in deaktivacijskih funkcij v odvisnosti od časa smo izvedli s pomočjo programskega orodja Mathematica ® za simbolično računanje. Ravnotežne potenciale 115 mV, 12 mV, 10.6 mV na K L E E E = =− = posameznih ionov smo izračunali z Nernstovo enačbo in v simulaciji uporabili maksimalne vrednosti prevodnosti 2 2 2 120 / , 36 / , 0.3 / Na K L G mS cm G mS cm G mS cm = = = ter vrednost membranske kapacitivnosti 2 1 / C F cm μ = . S tem lahko sistem enačb (8),(9) s predpisanimi začetnimi pogoji ( ) ( ) ( ) ( ) 0 , 0 , 0 in 0 V m h n zapišemo v obliki: ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 3 4 25 10 18 25 10 20 30 10 10 10 10 10 d 120 115 36 12 0.3 10.599 , d d 1 4 , d 1 d 0.07 1 , d 1 d 0.1 d 1 V t V t V t V t V t V t V t V t I t m h V t n V t V t t m m Exp m t Exp h h Exp h t Exp n t Exp − − − − − − − = − − − + − −                     = − −       −       = − −       +           =   −     ( ) ( ) 80 1 0.125 . V t n Exp n −                    − −         (12) in ga neposredno rešujemo z numerično integracijo. V simulaciji nas predvsem zanimata vpliv začetne depolarizacije membrane in vpliv stimulacijskega toka na širjenje akcijskega potenciala. Rezultati simulacije so prikazani na slikah 2, 3 in 4. Slika 2 prikazuje časovne poteke a) ene periode akcijskega potenciala (perioda T p =15 ms), b) vlaka impulzov akcijskega potenciala v daljšem časovnem razdobju T=100 ms, c) aktivacijskih funkcij m(t) in n(t) ter deaktivacijske funkcije h(t) in d) ionskih tokov I Na , I K in I L v času ene periode T p pri začetni depolarizaciji membrane V(0)=1 mV in stimulacijskem toku I=10 μA/cm 2 . Slika 3 prikazuje podobno zaporedje časovnih potekov akcijskega potenciala, vlaka impulzov akcijskega potenciala, aktivacijskih funkcij m(t) in n(t) ter deaktivacijske funkcije h(t) in ionskih tokov , in Na K L I I I pri povišani začetni depolarizaciji membrane V(0)=15 mV in stimulacijskem toku I=10 μA/cm 2 . Na sliki 4 pa je to zaporedje časovnih potekov prikazano pri začetni depolarizaciji membrane V(0)=1 mV in povišani vrednosti stimulacijskega toka I=50 μA/cm 2 . Rezultati kažejo, da večja depolarizacija membrane ne vpliva na trajanje periode akcijskega potenciala in maksimalno vrednost impulza (ki ohranjata vrednosti T p =15 ms, oziroma približno 100 mV), pač pa vpliva na njegovo zakasnitev, ki se zmanjša iz 6 ms na 2 ms. Razen tega se pri večji depolarizaciji membrane pojavi začetni padec akcijskega potenciala, ki mu sledi naraščajoča faza impulza. V padajoči fazi impulza se pojavi majhna grba, Anali PAZU - Letnik 1, leto 2011, številka 2 126 ki se prav tako časovno premakne (pri depolarizaciji z 1 mV se pojavi v času 7.5 ms, pri depolarizaciji s 15 mV pa v času 3.2 ms). Slika 2. Časovni poteki a) akcijskega potenciala, b) vlaka impulzov, c) vžignih spremenljivk kanalov in d) ionskih tokov kanalov pri stimulacijskem toku I=10 μA/cm 2 in začetni depolarizaciji membrane V(0)= 1 mV. Slika 3. Časovni poteki a) akcijskega potenciala, b) vlaka impulzov, c) vžignih spremenljivk kanalov in d) ionskih tokov kanalov pri stimulacijskem toku I=10 μA/cm 2 in začetni depolarizaciji membrane V(0)= 15 mV. 0 2 4 6 8 10 12 14 0 20 40 60 80 Čas t m s a Akcijski potencial V t  mV 0 20 40 60 80 100 0 20 40 60 80 Čas t m s b Akcijski potencial V t  mV 0 2 4 6 8 10 12 14 0.0 0.2 0.4 0.6 0.8 1.0 Čas t m s c Vžigne spremenljivke kanalov m t ,h t in n t n t h t m t 0 2 4 6 8 10 12 14  600  400  200 0 200 400 600 Čas t m s d Membranski tokovi ionov I Na  t , I K  t in I L  t I L  t I K  t I Na  t 0 2 4 6 8 10 12 14 0 20 40 60 80 100 Čas t m s a Akcijski potencial V t  mV 0 20 40 60 80 100 0 20 40 60 80 100 Čas t m s b Akcijski potencial V t  mV 0 2 4 6 8 10 12 14 0.0 0.2 0.4 0.6 0.8 1.0 Čas t m s c Vžigne spremenljivke kanalov m t ,h t in n t n t h t m t 0 2 4 6 8 10 12 14  600  400  200 0 200 400 600 800 Čas t m s d Membranski tokovi ionov I Na  t , I K  t in I L  t I L  t I K  t I Na  t Anali PAZU - Letnik 1, leto 2011, številka 2 127 Slika 4. Časovni poteki a) akcijskega potenciala, b) vlaka impulzov, c) vžignih spremenljivk kanalov in d) ionskih tokov kanalov pri stimulacijskem toku I=50 μA/cm 2 in začetni depolarizaciji membrane V(0)= 1 mV. Povečanje stimulacijskega toka vpliva na zmanjšanje periode T p akcijskega impulza, ki se iz 15 ms pri I=10 μA/cm 2 zniža na 9 ms pri I=50 μA/cm 2 . Značilno za povečanje stimulacijskega toka je tudi, da izzove povišano maksimalno vrednost pri prvem impulzu, nato pa se maksimalne vrednosti hitro zmanjšujejo proti stacionarni vrednosti. Aktivacijska funkcija m(t) se bistveno hitreje spreminja kot funkciji h(t) in n(t), zato funkcijo m(t) štejemo za hitro spremenljivko, funkciji h(t) in n(t) pa za počasni spremenljivki, kar je osnova FitzHugh-Nagumovega nizko-dimenzionalnega modela. Zanimive ugotovitve izhajajo tudi iz časovnih diagramov ionskih tokov, ki kažejo, da sta I Na in I K nasprotnega znaka, I L pa istega znaka kot I K , vendar je I L majhen v skladu s pričakovanji. 4. Zaključek V članku smo prikazali izpeljavo enačb Hodgkin- Huxleyevega modela celične membrane za vzdražljive (ekscitabilne) celice, kakršne so nevroni in mišične celice. V izpeljavi smo pokazali, da celično membrano lahko modeliramo kot električno vezje, pri čemer morajo imeti elementi vezja nekatere neobičajne lastnosti, ki posnemajo aktiviranje in deaktiviranje ionskih kanalov. Te lastnosti opišemo s kinetičnimi enačbami prvega reda. V predstavljeni simulaciji smo štiri enačbe Hodgkin- Huxleyevega modela reševali z numerično integracijo ob uporabi programskega orodja Mathematica. V simulaciji smo raziskali vpliv spreminjanja začetne depolarizacije membrane in povečevanja stimulacijskega toka na višino impulza, potek posameznih faz impulza in trajanje periode impulza. V nadalnjem delu upamo, da bomo z večstopenjsko homotopsko perturbacijsko metodo [6], ki smo jih razvili v dosedanjem raziskovalnem delu, raziskali globalne dinamične značilnosti širjenja akcijskega potenciala z uporabo nizko-dimenzionalnega modela. Reference 1. [1] Hodgkin, A. L., Huxley, A. F. “A quantitative description of membrane current and its application to conduction and excitation in nerve”. J. Physiol., 117: 500-544, 1952. 2. [2] FitzHugh, R. “Impulses and physiological states in theoretical models of nerve membrane”. Biophysical J. , 1: 445-466, 1961. 3. [3] Pušenjak, Rudi, Oblak, Maks, Tičar, Igor. “Nonstationary Vibration and Transition through Fundamental Resonance of Electromechanical Systems Forced by a Nonideal Energy Source“. International Journal for Nonlinear Sciences and Numerical Simulation, May 2009, vol. 10, no. 5, str. 637-660. 4. [4] Simitev, R. D., Biktashev, V. N. “Asymptotics of conduction velocity restitution in models of electrical excitation in the heart”. Bull Math Biol, 73(1): 72- 115, Jan 2011. 5. [5] Seydel, R. “Practical Bifurcation and Stability Anaysis. From Equilibrium to Chaos”. Second Edition. Springer-Verlag New York,Inc. 1994. 6. [6] Pušenjak, Rudi, Oblak, Maks. “Raziskave dinamike kaotičnih sistemov z uporabo večstopenjske homotopsko perturbacijske metode“. Slovensko društvo za mehaniko. Kuhljevi dnevi 2011, Zbornik del, Ljubljana , 203-210. 0 2 4 6 8 10 12 14 0 20 40 60 80 100 Čas t m s a Akcijski potencial V t  mV 0 20 40 60 80 100 0 20 40 60 80 100 Čas t m s b Akcijski potencial V t  mV 0 2 4 6 8 10 12 14 0.0 0.2 0.4 0.6 0.8 1.0 Čas t m s c Vžigne spremenljivke kanalov m t ,h t in n t n t h t m t 0 2 4 6 8 10 12 14  500 0 500 Čas t m s d Membranski tokovi ionov I Na  t , I K  t in I L  t I L  t I K  t I Na  t