The Hodgkin-Huxley Equation

Department of Biomedical Engineering,

Kyung Hee University

ABSTRACT

본 프로젝트에서는 세포막에서 시간의 흐름에 따라 변화하는 막전위와 [eqn: Na ^{+} ,`K ^{+}]등의 conductance의 변화를 시뮬레이션하기 위해서 신경세포를 전기적인 회로로 모델링한다. 이를 통해 전압의존성 [eqn: Na ^{+}]통로와 전압의존성 [eqn: K ^{+}]통로의 개방도에 관여하는 gate 변수 m,h,n을 설정하여 이에 대한 미분방정식과 세포막 전위에 대한 미분방정식을 얻을 수 있다. 이 두 가지 미분방정식 모두 numerical method 중에서 backward Euler method를 이용하여 [eqn: Na ^{+}]통로와 전압의존성 [eqn: K ^{+}]통로의 개방도에 관여하는 m(t), h(t), n(t)를 구하고 이를 통해 각각의 이온에 대한 conductance를 구하고 최종적으로 세포막 전위 v(t)에 대한 solution을 얻을 수 있었다. 이를 바탕으로 세포막 전위에 전기적 자극에 대한 조건을 바꿔가며 모델링한 세포막의 불응기와, 역치 전위를 확인할 수 있었고 hyperpolarizing pulse 조건에 따라 나타나는 결과를 통해 세포막에서 나타나는 ion conductance와 막전위의 관계에 대해 이해할 수 있었다.

TABLE OF CONTENTS

Abstract .......................................................................................................	1
Table Of Contents ........................................................................................	2
List of Figures and Tables ..........................................................................	3
Introduction .................................................................................................	4
Methods ....................................................................	5
Results and Discussion ................................................................................	13
Conclusions .................................................................................................	31
References ...................................................................................................	32
Appendices ..................................................................................................	33

LIST OF FIGURES AND TABLES

Figure 1	활동전위와 [eqn: Na ^{+} ,`K ^{+}] 전도도의 관계 그래프[3]	7
Figure 2	세포막의 Electrical Circuit Model [4]	7
Figure 3	전압이 step increase 후 step decrease할 때 [eqn: g _{K}]의 변화 [5]	9
Figure 4	전압이 A~J로 변화할 때 마다의 시간에 따른 [eqn: g _{Na}]의 변화 [5]	10
Figure 5	[eqn: alpha _{n}]의 값을 결정하는 alphan.m function의 flow chart	13
Figure 6	[eqn: beta _{n} `]의 값을 결정하는 betan.m function의 flow chart	14
Figure 7	[eqn: alpha _{m}]의 값을 결정하는 alpham.m function의 flow chart	15
Figure 8	[eqn: beta _{m} `]의 값을 결정하는 betam.m function의 flow chart	15
Figure 9	[eqn: alpha _{h}]의 값을 결정하는 alphah.m function의 flow chart	16
Figure 10	[eqn: beta _{m} `]의 값을 결정하는 betam.m function의 flow chart	16
Figure 11	snew function file의 flow chart	17
Figure 12	인가해준 전류의 변화	17
Figure 13	izero.m의 flow chart	18
Figure 14	voltage step increase	18
Figure 15	HH.m의 flow chart	20
Figure 16	활동전위의 발생과 10ms일 때 1ms동안 15uA의 전기적 자극에 의한 막전위, gate 변수 그래프	21
Figure 17	활동전위의 발생과 25ms일 때 1ms, 15uA의 전기적 자극에 의한 막전위, gate 변수 그래프	22
Figure 18	활동전위의 발생과 10ms일 때 2ms, 15uA의 전기적 자극에 의한 막전위, gate 변수 그래프	23
Figure 19	활동전위의 발생과 10ms일 때 1ms, 30uA의 전기적 자극에 의한 막전위, gate 변수 그래프	24
Figure 20	vstart=-70mV인 경우의 막전위와 gate variable그래프	25
Figure 21	vstart=-67mV 경우의 막전위와 gate variable그래프	25
Figure 22	vstart=-64mV일 경우의 막전위와 gate variable그래프	26
Figure 23	vstart=-63mV일 경우의 막전위와 gate variable그래프	26
Figure 24	vstart의 값의 변화(-73mV~-88mV)에 따른 막전위, gate variable, ion conductance의 그래프	28
Figure 25	vstart=-90mV일 경우의 막전위와 gate variable그래프	29
Figure 26	vstart=-90mV일 경우의 Na+와 K+의 conductance	29

Table 1	Hodgkin-Huxley 모델에 사용되는 변수들과 그 초기값의 의미	19

INTRODUCTION

신경세포란 신경계를 구성하는 세포이며 나트륨 통로, 칼륨 통로등의 이온 통로를 발현하여 다른 세포와는 달리 전기적인 방법으로 신호를 전달할 수 있다. 또한 인접한 다른 신경세포와 시냅스를 통해 화학적 신호를 주고받음으로써 다양한 정보를 받아들이고 저장하는 흥분성 세포이다. 본 프로젝트에서는 이러한 신경세포의 막의 구성, 물질운반, 전위, 흥분성 등의 특징으로 바탕으로 전기회로적으로 모델링하여 세포막에서 시간의 흐름에 따라 변화하는 막전위와 [eqn: Na ^{+} ,`K ^{+}]등의 conductance의 변화를 시뮬레이션 한다. 이러한 시뮬레이션은 1952년 Hodgkin–Huxley가 고안한 neuron에서의 활동전위를 묘사하는 수학적 모델인 Hodgkin–Huxley model을 기반으로 한다[Ref]. 이들은 squid giant axon에서의 활동전위가 발생하고 전파되는 것을 설명하기 위해 이 model을 고안하였고 이로 인해 1963년 노벨상을 수상하게 되었다[Ref]. Hodgkin–Huxley model에서는 세포막을 capcitor와 저항으로 구성한 전기회로로 단순화하여 1차 미분 방정식으로 나타낸다.

본 프로젝트에서는 세포막의 안팎을 이동하는 여러 가지 이온들 중에 [eqn: Na ^{+} ,`K ^{+}]의 이동에 따른 영향력이 크다는 것을 반영하여 이에 따르는 conductance를 계산하고 이들의 변화에 따른 막전위의 변화도 함께 살펴보았다. 이를 계산하기 위해서는 voltage clamp 방법을 사용하여 membrane potential을 고정하고 voltage constant를 유지하기 위해 흐르는 transient trans-membrane current를 측정한다. 이를 통해 시간의 흐름에 따른 막전위를 구하고 최종적으로 세포막을 이동하는 각각의 ion에 대한 conductance를 구할 수 있게 된다.

비록 생리학적으로 설명이 가능한 완벽한 모델을 구현한 것은 아니지만 몇 가지 가정과 실험적인 데이터를 바탕으로 fitting한 equation이 비교적 실제와 흡사하기 때문에 세포막에서 나타나는 활동전위에 대해서 이해하고 여러 가지 전기적 자극을 설정하여 시뮬레이션 해 봄으로써 세포막에서 일어나는 변화를 관찰할 수 있게 된다.

METHODS

본 프로젝트에서는 인체의 혈압변화를 알아보기 위한 시뮬레이션을 진행한다. 이때 systemic artery를 하나의 closed loop로 단순화하여 모델링하여 governing equation을 유도하였다.

Equipments

인체의 혈압변화 simulation을 위한 프로그램과 governing equation.

Conditions

1) 세포막의 구성 및 물질운반 기능

세포막은 인지질 이중층으로 구성되어 이온, 포도당, 요소와 같은 수용성 물질이 통과하는 것을 막고 산소, 이산화탄소, 알코올과 같은 지용성물질만을 통과시킨다. 또한 세포막에는 단백질로 된 pore들이 존재하여 여러 가지 물질들을 이동시키기거나 이들의 움직임을 조절한다. 세포막은 세포의 안과 밖을 구분하며 안과 밖 모두에 NaCl이나 KCl과 같은 염들이 물에 녹아 [eqn: Na ^{+} ,`K ^{+} ,`Cl ^{-}]과같은 이온들로 해리되어 존재한다. 이때 세포막은 세포의 안과 밖을 구분하는 경계가 되며 이 이온들을 선택적으로 이동시켜 세포 안팎의 농도를 일정하게 유지한다. Maybe add some references???

세포막을 통해서 여러 가지 물질들이 운반되는데 이는 세포막의 인지질층의 틈을 통해 이동하거나 단백질 pore를 통해 무작위적으로 이동하는 확산과 수송단백질을 통해 이온이나 다른물질들을 에너지 경사에 역행하여 이동시키는 능동수송에 의해 이루어진다. 확산에는 단순확산과 촉진확산이 있으며 이중에서 단순확산은 지용성물질을 인지지질 이중막의 틈새로 이동시키는 경우와 [eqn: Na ^{+} channel과``K ^{+} `channel]과 같은 protein channel을 통해 선택적으로 이동시키는 경우가 있으며 촉진확산은 carrier protein과의 결합을 통해 glucose와 같은 물질들을 이동시킨다. 능동수송에는 [eqn: Na ^{+} -K ^{+} `pump]를 통해 소듐이온을 세포밖으로 내보내는 동시에 포타슘 이온을 세포 안으로 옮기는 것이 있으며 그 외에도 칼슘이온, 수소이온 등을 pump를 통해 농도 경사에 역행하여 이동시킨다[1]. 본 프로젝트에서는 이와 같은 여러 가지 물질이동을 통해 생성된 세포막의 전위와 활동전위를 시뮬레이션한다.

2) 세포막 전위

앞서 설명하였듯 세포 내외에서 [eqn: Na ^{+} ,`K ^{+} ,`Cl ^{-}]과 같은 이온들이 확산과 능동수송에  의해 이동하면서 전위가 발생하게 된다. 평형상태에서의 각각의 이온의 확산에 의한 세포막 전위차는 다음과 같은 Nernst Potential로 나타낼 수 있다.

[eqn: V _{s} `=` {RT} over {ZF} ln( {[S] _{o}} over {[S] _{i}} )]                                 (1)

(여기서, R은 기체 상수, T는 절대 온도, F는 Faraday 상수, Z는 이온S의 전하량, [eqn: [s] _{o}]은 세포막 밖 이온s의 농도, [eqn: [s] _{i}]은 세포막 안 이온 S의 농도이다.)

이를 각각의 이온이 가지는 극성, 막투과성, 농도를 반영하여 합산하는 방식이 Goldman-Hodgkin-Katz Equation이며 다음과 같이 나타난다.

[eqn: V _{r} `=- {RT} over {F} ln( {P _{Na} [Na ^{+} ] _{i} +P _{K} [K ^{+} ] _{i} +P _{Cl} [Cl ^{-} ] _{o}} over {P _{Na} [Na ^{+} ] _{o} +P _{K} [K ^{+} ] _{o} +P _{Cl} [Cl ^{-} ] _{i}} )]                (2)

(여기서 R은 기체 상수, T는 절대 온도, F는 Faraday 상수, [eqn: P _{ion}]는 각이온에 대한 투과성, [eqn: [ion] _{i}]은 세포 내부의 이온 농도, [eqn: [ion] _{o}]은 세포 외부의 이온 농도를 의미한다.) 본 프로젝트에서는 이와 같이 [eqn: Na ^{+} ,`K ^{+} ,`Cl ^{-}]이온에 의해 전위차가 발생하는 것을 이해하고 이 식을 이용하여 세포막의 전위와 활동전위를 시뮬레이션한다.

3) 세포의 흥분성

세포는 ionic channel을 조절하여 막전위를 조절하는 아주 중요한 역할을 한다. 특히 신경세포와 근육세포는 membrane potential을 신호로 하여 이를 발생시키고 전달하여 신경계의 작동이나 근육의 수축을 일으킨다. 세포에는 내장의 표피세포나 photoreceptor과 같이 막전위는 발생하지만 활동전위가 발생하지 않는 비 흥분성 세포들과 심장, 근육, 신경세포와 같이 활동전위를 발생시키는 흥분성 세포가 있다. 본 프로젝트에서 선택한 모델은 바로 흥분성 세포인 신경세포이다[1].

근육, 신경 등 흥분성 세포의 흥분에 의한 세포막의 일시적인 전위변화를 ‘활동전위(action potential)’이라고 한다. 세포는 [eqn: Na ^{+} ,`K ^{+}]등의 여러 이온 펌프의 활동에 의해 세포 안팎의 이온 조성 차이로 –90mV~-60mV의 안정막 전위를 가진다.

한편 흥분성 세포가 흥분을 하게 되면 세포막 전위가 –70mV~-50mV가 되면서 전압의존성 [eqn: Na ^{+}]통로를 활성화(m-gate)시켜 [eqn: Na ^{+}]통로를 통해 [eqn: Na ^{+}]가 세포 안으로 유입된다. 이때 세포막 내부는 +30mV~+40mV의 전위를 가지게 되며 이 현상을 탈분극(depolarization)이라고 한다.

전압의존성 [eqn: Na ^{+}]통로는 활성화(m-gate)됨과 동시에 불활성화(h-gate)가 일어나며 세포막 전위가 –90mV~0mV에 이르기까지 전압의존성 [eqn: K ^{+}]통로(n-gate)의 문이 서서히 열린다. 즉, [eqn: Na ^{+}]통로의 불활성화 될 때 ([eqn: Na ^{+}]이온의 유입이 감소할 때) [eqn: K ^{+}]통로(n-gate)가 활성화 되어 [eqn: K ^{+}]가 세포밖으로 빠져나가게 된다. 이를 통해 세포막 전위는 다시–90mV~-60mV의 안정막 전위를 회복하게 되는데 이를 재분극(repolarization)이라고 한다.

이러한 재분극 과정에서 [eqn: K ^{+}]통로(n-gate)가 서서히 닫히기 때문에 포타슘 이온이 세포 밖으로 많은 양이 빠져나가서 본래 안정막 전위보다 더 낮은 전위(–70mV~–75 mV)를 가지는데 이를 과분극(hyperpolarization)이라고 한다[2].

이처럼 소듐과 포타슘 통로의 열리고 닫힘에 의해서 활동전위는 발생하게 되며 이를 그래프로 나타낸 것이 Figure 1이다. 본 프로젝트에서는 이와 같이 [eqn: Na ^{+} ,`K ^{+}]등의 여러 이온들의 conductance와 활동전위의 관계를 이해하고 이를 시뮬레이션한다.

Figure 1.	활동전위와 [eqn: Na ^{+} ,`K ^{+}] 전도도의 관계 그래프[3]

3. Governing equation

1) 세포막의 Electrical Circuit Model

세포막은 전하를 세포 내외로 분리시키기 때문에 Capacitor로 나타낼 수 있으며 인지질 이중층으로 구성되어 물질을 선택적으로 투과하기 때문에 이 성질을 저항으로 나타낼 수 있다. 따라서 이를 전기적 회로 모델로 나타내면 다음 그림과 같다.

Figure 2.	세포막의 Electrical Circuit Model [4]

커패시터의 특성을 반영하여 [eqn: C _{m}]에 관계식은 다음과 같이 세울 수 있다.

[eqn: C _{m} `=` {Q} over {V _{c}}], [eqn: C _{m} V _{c} `=`Q], [eqn: C` _{m} {dv} over {dt} `=` {dQ} over {dt} `=`I _{c}]                  (3)

이를 키르히호프 법칙에 적용하여 식을 세우면 다음과 같이 나타낼 수 있다.

[eqn: C _{m} {dv} over {dt} `+`I _{ion} `=`0]                                (4)

이 회로에서 농도 차에 의해 발생하는 전위차는 앞서 정리했던 식(1)이고, 전류에 의해 발생하는 전위차는 다음과 같이 나타낼 수 있다.

[eqn: rI _{ion}]                                       (5)

(여기서, r은 세포막 통로의 저항, [eqn: I _{ion}]는 세포막 안팎을 통과하는 전류이다.)

따라서 이 두 가지 전위차를 더하여 세포막을 안팎을 통과하는 전류 [eqn: I _{ion}]에 대해 정리하면 다음과 같다.

[eqn: v`=`rI _{ion} `+`E _{s}]                                  (6)

[eqn: I _{ion} `=` {1} over {r} (v-E _{s} )`=`g(v-E _{s} )]                            (7)

(여기서 g는 세포막의 conductance이다.)

따라서 최종적으로 식(4)에 식(7)을 대입하여 얻어진 식은 다음과 같다.

[eqn: C _{m} {dv} over {dt} `+`g _{s} (v`-`E _{s} )`=`0]                             (8)

2) Hodgkin-Huxley Model

본 프로젝트에서 simulation할 것은 신경세포이기 때문에 주된 ionic 전류는 [eqn: Na ^{+} ,`K ^{+}]에 의한 것으로 그 외에는 [eqn: Cl ^{-}]로 대표되는 여러 이온들에 의한 아주 작은 leakage current가 있다. 이를 식(8)에 적용하면 다음과 같은 미분방정식이 나타난다.

[eqn: C _{m} {dv} over {dt} `=`-g _{Na} (v-E _{Na} )`-`g _{K} (v-E _{K} )`-`g _{C`l} (v-E _{C`l} )`+`I _{app}]            (9)

(여기서, [eqn: I _{app}]는 인가해준 전류이다.)

이를 단순화하면 다음과 같은 간단한 미분방정식을 얻을 수 있다.

[eqn: C _{m} {dv} over {dt} `=`-g _{eff} (v-E)`+`I _{app}]                       (10)

(여기서 [eqn: g _{eff} `=`g _{Na} `+`g _{K} `+`g _{C`l}], [eqn: E=(g _{Na} E _{Na} +g _{K} E _{K} +g _{L} E _{L} )/g _{eff}]로 안정막전위를 의미한다.)

한편 인가한 전류가 steady current라면, 전압의 값이 시간에 따라 변하지 않기 때문에 [eqn: {dv} over {dt}]=0가 된다. 이를 반영하여 전압 V에 대한 식을 정리하면 다음과 같다.

[eqn: v`=`E`+ {I _{app}} over {g _{eff}}]                             (11)

3) 시간에 따라 변화하는 Conductance - Voltage clamp방법

충분히 작은 전류를 인가해준다면 위의 식 (11)은 실제 현상을 잘 설명한다고 볼 수 있지만 큰 전류가 인가된다면 conductance가 시간의 흐름에 따라 일정하지 않기 때문에 오차가 생기게 된다. Hodgkin-Huxley는 이를 해결하기 위해 voltage clamp 방법을 사용하였다. voltage clamp란 membrane potential을 고정하고 voltage constant를 유지하기 위해 흐르는 전류를 측정한다. 이때 측정된 전류는 transient transmembrane current이다.

이를 통해 몇 가지 무척추 동물의 큰 신경 섬유를 불리하여 실험을 진행할 수 있게 되었는데 이때 세포막 바깥에 tetrodotoxin(복어독)을 적용하여 sodium channel을 억제하거나 세포막 안에 tetraethylammonium ion을 적용하여 potassium channel을 억제하는 등의 방법을 통해 각각의 channel을 통한 전류를 측정할 수 있었다. 그 결과 처음에 발생한 안쪽으로 향하는 전류는 [eqn: Na ^{+}]에 의한 것이고 나중에 나타나는 밖으로 향하는 전류는 [eqn: K ^{+}]에 의한 것임을 알 수 있었다. 이 중에서 먼저 [eqn: K ^{+}]의 conductance부터 살펴보도록 하자.

① Potassium Conductance

다음의 Figure 3은 실험적인 결과를 통해 전압이 step increase 후 step decrease할 때 [eqn: g _{K}]의 변화를 나타낸 그래프이다. 이 그래프의 앞부분 (a)에서는 sigmoidal increase가 나타나고 그 이후(b)에서 exponential decrease가 나타났다.

Figure 3.	전압이 step increase 후 step decrease할 때 [eqn: g _{K}]의 변화 [5]

Hodgkin과 Huxely는 이러한 실험적인 데이터를 통해 conductance에 대한 식을 나타내고자 했다. 하지만 이를 생리학적으로 설명할수 없어 그래프의 fitting을 위하여 [eqn: K ^{+}]통로는 4개의 n-gate를 가지고 각각의 gate는 동일하며 독립적으로 열리거나 닫힐 수 있다는 가정을 하였고 다음과 같은 관계식으로 potassium conductance를 정의하였다.

[eqn: g _{K} `=` {bar{g _{K} `}} `n ^{4}]                                (12)

(여기서, [eqn: {bar{g _{K}}}]는 상수(=36)이고, n은 [eqn: K ^{+}]통로의 gate 활성도이다.)

즉, potassium conductance는 변수 n의 네 제곱으로 표현이 가능하며 비례관계를 가진다는 것을 알 수 있다. 여기서 변수 n이 gate가 열린 상태일 확률을 나타내고 [eqn: alpha _{n}]은 1초 동안 닫힌 상태의 gate가 열리는 횟수, [eqn: beta _{n} `]은 1초동안 열린 상태의 gate가 닫히는 횟수를 의미한다고 할 때 다음과 같은 n에 관한 미분 방정식을 세울 수 있다.

[eqn: {dn} over {dt} `=` alpha _{n} (1-n)`-` beta _{n} `n]                            (15)

한편 potassium conductance는 Figure 3에서 나타난 바와 같이 증가하는 부분(a)와 감소하는 부분(b)가 나타나는데 이에 맞게 fitting하여 해당하는 식을 찾을 수 있다. 먼저 sigmoidal increase가 나타나는 Figure2의 (a)부분에 맞게 미분방정식을 풀면

[eqn: n(t)=``n _{INF } (v)[1-exp( {-t} over {tau _{n} (v _{o} )} )]],                      (16)

한편 exponential decrease가 나타난 Figure2의 (b)부분에 맞게 미분방정식을 풀면

[eqn: n(t)=``n _{INF } (v)exp( {-t} over {tau _{n} (v _{o} )} )].                         (17)

여기서 [eqn: n _{INF } (v)= {alpha _{n} (v)} over {alpha _{n} (v)+ beta _{n} (v)}]이고 [eqn: tau _{n} (v)= {1} over {alpha _{n} (v)+ beta _{n} (v)}]을 의미하기 때문에 이를 이용하여 [eqn: alpha _{n} `,` beta _{n}]을 구하면 다음과 같다.

[eqn: alpha _{n} =0.1 {(v+60)/10} over {1-exp(-(v+60)/10)}]                       (18)

[eqn: beta _{n} =0.125exp(- {v+70} over {80} )]                          (19)

② Sodium Conductance

한편 sodium conductance는 앞서 potassium conductance를 구했던 방법과 동일하게 실험적으로 얻은 그래프를 통해 식을 세울 수 있지만 더 많은 변수가 필요하다. 우선 실험적으로 나타난 sodium conductance의 그래프는 다음과 같다.

Figure 4.	전압이 A~J로 변화할 때 마다의 시간에 따른 [eqn: g _{Na}]의 변화 [5]

sodium conductance의 그래프도 생리학적으로 설명할 수 없기 때문에 Hodgkin-Huxley는 potassium conductance를 구하는 방법과 마찬가지로 sodium channel이 4개의 gate를 가지고 있으며 이 4개의 gate들은 potassium channel과 달리 모든 gate가 동일하지 않고 3개의 m-gate와 1개의 h-gate를 가지고 있다고 가정하였다. 이들은 독립적으로 열리거나 닫힐수 있으며 m-gate는 전압이 증가갈수록 channel이 열리는 데에 기여하지만 h-gate는 전압이 증가할수록 channel이 닫히는 데에 기여한다. 따라서 sodium conductance는 다음과 같은 관계식으로 나타낼 수 있다.

[eqn: g _{Na} (v)`=` {bar{g _{Na}}} m ^{3} h]                               (20)

(여기서, [eqn: {bar{g _{Na}}}]는 상수(=120)이고, m은 [eqn: Na ^{+}]통로의 gate 활성도이고, h는 [eqn: Na ^{+}]통로의 gate 비활성도이다.)

potassium conductance를 구했던 방법과 마찬가지로 변수 m과 h가 각각의 gate가 열린 상태일 확률을 나타내고 sodium channel의 m-gate와 h-gate도 1초 동안 닫힌 상태의 gate가 열리는 횟수 [eqn: alpha _{m}]과 [eqn: alpha _{h}], 1초 동안 열린 상태의 gate가 닫히는 횟수 [eqn: beta _{m} `]과 [eqn: beta _{h} `]를 이용하여 다음과 같은 m과 h에 관한 미분방정식을 세울 수 있다.

[eqn: {dm} over {dt} `=` alpha _{m} (1-m)- beta _{m} m]                            (21)

[eqn: {dh} over {dt} `=` alpha _{h} (1-h)`-` beta _{h} `h]                             (22)

이 또한 실험적으로 얻어진 Figure 4와 같은 그래프에 line fitting하면 최종적으로 다음과 같이 [eqn: alpha _{m}], [eqn: alpha _{h}], [eqn: beta _{m} `], [eqn: beta _{h} `]에 대한 식을 얻을 수 있다.

[eqn: alpha _{m} =1.0 {(v+45)/10} over {1-exp(-(v+45)/10)}]                      (23)

[eqn: beta _{m} =4.0`exp(-(v+70)/18)]                         (24)

[eqn: alpha _{h} =0.07`exp(-(v+70)/20)]                        (25)

[eqn: beta _{h} =1.0 {1} over {1+exp(-(v+60)/10)}]                       (26)

③ Leakage conductance

Hodgkin-Huxley는 [eqn: Na ^{+} ,`K ^{+}]channel 외의 [eqn: Cl ^{-}]로 대표되는 여러 이온의 channel을 찾아냈고 이를 leakage channel이라고 한다. [eqn: Na ^{+} ,`K ^{+}]channel과 달리 leakage channel의 conductance는 다음과 같이 항상 일정한 값의 상수를 나타냈다.

[eqn: g _{L} `=` {bar{g _{L}}}]                                   (27)

4) 최종적인 governing equation

앞서 [eqn: Na ^{+} ,`K ^{+}]의 coductance와 관련하는 변수 m,h,n에 대한 미분방정식을 얻어냈다면 최종적으로 numerical method를 통해 solution을 구하는 과정이 필요하다. 다음의 식 (28)에서 변수 s를 m,h,n을 대표하는 변수로 설정한다면 numerical solution 중에서도 Backward Eurler Method를 사용하여 다음과 같이 나타낼 수 있다.

[eqn: {s(t)-s(t- TRIANGLE t)} over {TRIANGLE t} = alpha _{s} (v(t- TRIANGLE t))(1-s(t))- beta _{s} (v(t- TRIANGLE t))s(t)]        (28)

이 미분방정식의 해를 구하면 다음과 같은 함수를 얻을 수 있다.

[eqn: s(t)= {s(t- TRIANGLE t)+ TRIANGLE t alpha _{s} v(t- TRIANGLE t)} over {1+ TRIANGLE t( alpha _{s} v(t- TRIANGLE t)+ beta _{s} v(t- TRIANGLE t))}]                  (29)

여기서 전압이 v가 일정하게 유지된다고 했을 때 변수 s는 다음과 같은 관계식 (30)으로 정의될 수 있다.

[eqn: s= {alpha _{s} (v)} over {alpha _{s} (v)+ beta _{s} (v)}]                                  (30)

또한 앞서 구했던 막전위에 관한 미분방정식 식(10) 또한 Backward Eurler Method를 사용하여 다음과 같이 나타낼 수 있다.

[eqn: C _{m} {v(t)-v(t- TRIANGLE t)} over {TRIANGLE t} +g _{eff} (t)(v(t)-E(t))=I _{app} (t)]               (31)

이 미분방정식을 해를 구하면 다음과 같은 함수를 얻을 수 있다.

[eqn: v(t)=` {v(t)`+` {TRIANGLE t} over {C} (g(t)E(t)+`I _{app} (t))} over {1+` {TRIANGLE t} over {C} g(t)}]                 (32)

이때 다음과 같은 parameter constant 값들을 대입하여 최종적으로 세포막 전위에 대한 시뮬레이션을 진행한다.

[eqn: C _{m} =1.0`[ mu A TIMES ms/cm ^{2} ]]                           (33)

[eqn: {bar{g}} _{Na} =120`[ mu A/(mV BULLET cm ^{2} )]], [eqn: {bar{g}} _{K} =36`[ mu A/(mV BULLET cm ^{2} )]], [eqn: {bar{g}} _{L} =0.3`[ mu A/(mV BULLET cm ^{2} )]] (34)

[eqn: E _{Na} =45`[mV]], [eqn: E _{K} =-82`[mV]], [eqn: E _{L} =-59`[mV]]                (35)

RESULTS AND DISCUSSIONS

(Exercise 1) Hodgkin-Huxley equation을 시뮬레이션 하기 위해 사용한 file을 분석해보자.

1) alphan.m

앞서 governing equation에서 설명하였듯 [eqn: alpha _{n}]은 1초 동안 닫힌 상태의 n-gate([eqn: `K ^{+}]통로의 gate)가 열리는 횟수를 의미하며 line fitting을 통해 얻은 식이 바로 식 (18)이다.

[eqn: alpha _{n} =0.1 {(v+60)/10} over {1-exp(-(v+60)/10)}]                       (18)

이 식에서 값이 정의되지 않는 부분은 분모가 0이 되는 순간이며 이는 [eqn: {v+60} over {10}]가 0이 될 때이다. 이때의 [eqn: alpha _{n}] 값을 로피탈의 법칙에 따라 0.1이 되도록 설정하고 이외의 경우에 식 (18)에 따라 [eqn: alpha _{n}]를 v값을 입력받아 계산하도록 하여 변수 a에 저장하는 function파일이 alphan,m이다. 편의를 위해 [eqn: {v+60} over {10}]를 theta로 정의하였고 이에 따른 function file의 flow chart는 Figure 5와 같다.

Figure 5.	[eqn: alpha _{n}]의 값을 결정하는 alphan.m function의 flow chart

2) betan.m

한편 [eqn: beta _{n} `]은 1초 동안 열린 상태의 n-gate([eqn: `K ^{+}]통로의 gate)가 닫히는 횟수를 의미하며 line fitting을 통해 얻은 식이 바로 식 (19)이다.

[eqn: beta _{n} =0.125exp(- {v+70} over {80} )]                          (19)

모든 범위의 v에 의해 [eqn: beta _{n} `]이 결정되어 그 값을 변수 b에 저장하는 function file이다. 편의를 위해 [eqn: {v+70} over {80}]를 theta로 정의하였고 이를 반영한 betan.m의 flow chart는 Figure 6와 같다.

Figure 6.	[eqn: beta _{n} `]의 값을 결정하는 betan.m function의 flow chart

3) alpham.m

[eqn: alpha _{m}]은 1초 동안 닫힌 상태의 m-gate([eqn: `Na ^{+}]통로의 activation gate)가 열리는 횟수를 의미하며 line fitting을 통해 얻은 식이 바로 식 (23)이다.

[eqn: alpha _{m} =1.0 {(v+45)/10} over {1-exp(-(v+45)/10)}]                      (23)

이 식에서 값이 정의되지 않는 부분은 분모가 0이 되는 순간이며 이는 [eqn: {v+45} over {10}]가 0이 될 때이다. 이때의 [eqn: alpha _{m}] 값을 로피탈의 법칙에 따라 1.0이 되도록 설정하고 이외의 경우에 식 (18)에 따라 [eqn: alpha _{m}]를 v값을 입력받아 계산하도록 하여 변수 a에 저장하는 function파일이 alpham,m이다. 편의를 위해 [eqn: {v+45} over {10}]를 theta로 정의하였고 이에 따른 function file의 flow chart는 Figure 7과 같다.

Figure 7.	[eqn: alpha _{m}]의 값을 결정하는 alpham.m function의 flow chart

4) betam.m

[eqn: beta _{m} `]은 1초 동안 열린 상태의 m-gate([eqn: `Na ^{+}]통로의 activation gate)가 닫히는 횟수를 의미하며 line fitting을 통해 얻은 식이 바로 식 (24)이다.

[eqn: beta _{m} =4.0`exp(-(v+70)/18)]                         (24)

모든 범위의 v에 의해 [eqn: beta _{m} `]이 결정되어 그 값을 변수 b에 저장하는 function file이다. 편의를 위해 [eqn: {v+70} over {18}]를 theta로 정의하였고 이를 반영한 betam.m의 flow chart는 Figure 8 같다.

Figure 8.	[eqn: beta _{m} `]의 값을 결정하는 betam.m function의 flow chart

5) alphah.m

[eqn: alpha _{h}]은 1초 동안 닫힌 상태의 h-gate([eqn: `Na ^{+}]통로의 inactivation gate)가 열리는 횟수를 의미하며 line fitting을 통해 얻은 식이 바로 식 (25)이다.

[eqn: alpha _{h} =0.07`exp(-(v+70)/20)]                        (25)

모든 범위의 v에 의해 [eqn: alpha _{h}]가 결정되어 그 값을 변수 a에 저장하는 function file이다. 편의를 위해 [eqn: {v+70} over {20}]를 theta로 정의하였고 이에 따른 function file의 flow chart는  같다.

Figure 9.	[eqn: alpha _{h}]의 값을 결정하는 alphah.m function의 flow chart

6) betah.m

[eqn: beta _{h}]은 1초 동안 닫힌 상태의 h-gate([eqn: `Na ^{+}]통로의 inactivation gate)가 닫히는 횟수를 의미하며 line fitting을 통해 얻은 식이 바로 식 (26)이다.

[eqn: beta _{h} =1.0 {1} over {1+exp(-(v+60)/10)}]                       (26)

모든 범위의 v에 의해 [eqn: beta _{m} `]이 결정되어 그 값을 변수 b에 저장하는 function file이다. 편의를 위해 [eqn: {v+60} over {10}]를 theta로 정의하였고 이를 반영한 betah.m의 flow chart는 같다.

Figure 10.	[eqn: beta _{m} `]의 값을 결정하는 betam.m function의 flow chart

7) snew.m

이 파일에서는 [eqn: Na ^{+}]와 [eqn: K ^{+}]channel의 conductance와 관련하는 변수인 n,m,h를 값을 새로운 값으로 update하는 과정이 일어난다. m의 경우를 예를 들어 설명해보면 m의 이전 값이 s_old로 전해지고 앞서 3)~4)의 function file을 통해 최종적으로 얻은 a,b가 alpha와 beta로 전해진다. 마지막으로 dt가 in_HH.m file로부터 전해짐으로써 입력을 받고 이를 통해 새로운 값 s를 앞서 구했던 governing equation 식 (29)를 이용하여 구할 수 있다.

[eqn: s(t)= {s(t- TRIANGLE t)+ TRIANGLE t alpha _{s} v(t- TRIANGLE t)} over {1+ TRIANGLE t( alpha _{s} v(t- TRIANGLE t)+ beta _{s} v(t- TRIANGLE t))}]                  (29)

마지막으로 이 file은 self-checking 과정을 포함한다. global 변수 check의 값이 1일 경우 함수가 snew가 호출될때마다 chsnew값이 출력되며 특별한 경우가 아닌 이상 [eqn: +- 0.5]범위 이내에 있어야 한다. Figure 11는 이와 같은 snew function file의 과정을 flow chart로 나타낸 것이다.

Figure 11.	snew function file의 flow chart

8) izero.m

본 시뮬레이션에서는 신경세포에서의 action potential을 발생시키기 위해 t1p와 t2p사이의 시간 동안 ip라는 전류가 흐르는 구형파의 전류를 인가한다(Figure 12). 이러한 전류의 변화를 나타낸 것이 바로 izero.m file이다. Figure 13는 이를 flow chart로 나타낸 것이다.

Figure 12.	인가해준 전류의 변화

Figure 13.	izero.m의 flow chart

9) in_mhnv.m

in_mhnv.m 파일은 m,h,n의 값과 v가 가지는 초기값을 설정해준다. 이 과정에서 우리는 전류를 인가함에 따라 전압이 변화하는 순간을 t=0으로하여 이 시점 이전에 안정막 전위인 –70mV를 holding voltage, vhold로 설정해준다. 이에 따라 전류가 인가되면 t=0인 순간에 전압이 증가하여 특정한 전압 vstart이 될 것이다. 하지만 이 충격은 아주 잠시동안 지속되기 때문에 Gate 변수 m,h,n은 변화할 시간이 없게 된다. 따라서 우리는 Figure 14에서 나타난 바와 같이 전압이 증가하도록 설정하고 이에 따라 v=vhold인 구간에서의 변하지 않는 m,n,h의 초기값을 앞서 정의한 governing equation (30)를 이용하여 구할 수 있다. 이를 통해 m,n,h의 초기값을 설정한 다음 최종적으로 v값을 vstart로 만들어주는 것이 이 파일의 역할이다.

Figure 14.	voltage step increase

10) in_HH.m

이 파일은 Hodgkin-Huxley 모델에 사용되는 변수들을 초기화 시키기 위한 값들을 지정하는 파일이다. 각각의 변수들의 이름의 의미와 초기값은 Table 4에 나타나 있다.

Tabel 1	Hodgkin-Huxley 모델에 사용되는 변수들과 그 초기값의 의미
	Variable	초기값	[Unit]	의미
	C	1.0	[eqn: LEFT [ mu A TIMES ms/cm ^{2} RIGHT ]]	세포막을 의미하는 커패시터의 용량값
	gNabar ([eqn: {bar{g _{Na}}}])	120	[[eqn: ( mu A/mV)/cm ^{2}]]	[eqn: Na ^{+}]통로의 모든 gate가 열렸을 때의 [eqn: Na ^{+}]conductance
	gKbar ([eqn: {bar{g _{K}}}])	36	[[eqn: ( mu A/mV)/cm ^{2}]]	[eqn: K ^{+}]통로의 모든 gate가 열렸을 때의 [eqn: K ^{+}]conductance
	gLbar ([eqn: {bar{g _{L}}}])	0.3	[[eqn: ( mu A/mV)/cm ^{2}]]	leakage channel의 conductance
	ENa ([eqn: E _{Na}])	45	[mV]	[eqn: Na ^{+}] channel의 평형 전위
	EK ([eqn: E _{K}])	-82	[mV]	[eqn: K ^{+}] channel의 평형 전위
	EL ([eqn: E _{L}])	-59	[mV]	leakage channel의 평형 전위
	dt	0.1	[ms]	시간 간격
	tmax	50	[ms]	시뮬레이션 진행 시간
	klokmax	tmax/dt	[개수]	시간 간격의 총 개수
	mhn_plot	zeros(3, klockmax)	[[eqn: 3` TIMES `klockmax]]	m,h,n 데이터(3가지)를 저장할 배열 생성
	v_plot	zeros(1, klockmax)	[[eqn: 3` TIMES `klockmax]]	v에 대한 데이터를 저장할 배열 생성
	t_plot	zeros(1, klockmax)	[[eqn: 3` TIMES `klockmax]]	t에 대한 데이터를 저장할 배열 생성
	vhold	-70	[mV]	신경 세포막의 전위를 –70[mV]로 고정(t = 0)
	vstart	-55	[mV]	전류 인가 후의 전위 (t > 0) -55mV는 활동전위를 일으키기에 충분
	t1p	10	[ms]	전류 인가 시작 시간
	t2p	11	[ms]	전류 인가 종료 시간
	ip	15	[[eqn: mu A]]	인가한 currentd의 크기
	check	1	[no unit]	check의 값이 1일 때 self-check 진행하여 chsnew 출력

11) HH.m

이 파일은 먼저 snew.m file을 이용하여 m(t), n(t), h(t)를 새로운 값으로 갱신하고 이에 따라 결정되는 [eqn: Na ^{+}]conductance([eqn: g _{Na}])와 [eqn: K ^{+}]conductance([eqn: g _{K}])를 governing equation (12),(20)을 통해 결정하고 izero.m file에서 발생시킨 전류를 반영하여 최종적으로 V에 대한 governing equation (32)을 이용하여 v값을 갱신한다. 그다음 self-check 과정을 거처 시간에 따른 m,n,h 데이터를 저장하고 이를 ploting한다. 이 파일의 flow chart는 Figure 15와 같다.

Figure 15.	HH.m의 flow chart

(Exercise 2) 신경 세포막에서의 활동전위 발생

HH.m를 실행함에 따라 세포막에는 처음에 –70mV의 전압(vhold)이 주어졌고 t=0이 되는 순간에 –55mV의 전압(vstart)이 될 수 있도록 전압을 만들어주었다. 이에 따라 세포막에서는 활동전위가 발생하게 되고 그에 따른 그래프가 0ms~10ms 사이의 그래프이다. 우선 세포막의 전위가 –70mV~-50mV가 되면서 [eqn: Na ^{+}]통로의 활성화와 관련된 m의 값은 0에 근접한 값에서 0~2.3ms안에 1까지 아주 급격히 증가하고 2.3~5ms사이에 급격히 감소하는 양상을 보였다. 한편 [eqn: Na ^{+}]통로의 불활성화와 관련된 h의 값은 0.6에서 시작하여 3ms일 때 최소값 0.07701이되고 다시 증가하여 10ms에는 0.5의 값을 가지게 되었다. [eqn: K ^{+}]통로와 관련된 n은 0~3ms사이에서는 0.1에서 0.7736까지 급격하게 증가하는 편이지만 그 이후에는 완만하게 감소하여 10ms에서는 0.3741의 값을 가진다. 한편 10ms와 11ms 사이에 전류 ip[mA]를 새롭게 인가하였기 때문에 이 부근에서 membrane potential을 나타내는 그래프에서 –76mV에서 –66mV 정도로 미세하게 증가하는 것을 확인할 수 있었다. 이에 따라 m은 약 0.02에서 0.06으로, n은 약0.37에서 0.365으로 ,h는 약0.51~0.52로 아주 미세한 차이가 발생했을 뿐 큰 차이가 발생하지 않았다.

이러한 결과는 세포막의 전위가 –70mV에서 –50mV로 증가할 때 전압의존성 [eqn: Na ^{+}]통로의 활성화가 진행되기 때문에 이와 관련한 m-gate가 열리게 되어 그래프에서 나타난 바와 같이 급격하게 m-gate의 개방도를 의미하는 m의 값이 증가하게 된다. 이에 따라 [eqn: Na ^{+}]는 세포 내로 유입되어 membrane potential을 높이는데 기여하여 막전위를 약 30mV까지 높인다. 전압의존성 [eqn: Na ^{+}]통로는 빠르게 반응이 일어나는 통로로써 m의 값은 급격히 감소하게 되고 활성화와 동시에 불활성화문이 닫히기 때문에 이와 관련된 h(h-gate의 개방도를 나타냄)의 값이 감소하게 된다(이 또한 빠르게 다시 원래의 상태를 회복하여 h값은 다시 증가). 한편 전압의존성 [eqn: K ^{+}]통로는 n-gate만이 관여하며 전압의존성 [eqn: Na ^{+}]통로와 달리 서서히 열리고 서서히 닫히는 특징이 있다. 이는 그래프에서 나타난 현상과 일치하는 것이며 안정 상태에서 닫혀 있던 전압의존성 [eqn: K ^{+}]통로가 막전위가 –90mV에서 0mV 사이에서 활성화되어 서서히 [eqn: K ^{+}]통로가 열려 [eqn: K ^{+}]가 세포밖으로 유출된다. [eqn: K ^{+}]통로는 다시 서서히 닫히면서 세포막 전위를 처음의 안정막전위 –70mV로 서서히 회복시킨다. 따라서 이러한 이론에 따라 0~10ms사이의 그래프는 생리학적으로 적합한 시뮬레이션이라고 볼 수 있다. 10ms에 1ms동안 15[eqn: mu A]의 전류를 인가하였기 때문에 약간의 전위가 상승하는 모습을 보이지만 앞서 나타난 것과 같은 활동전위가 발생하지 않은 이유는 전기 자극이 아주 약하기 때문으로 추정된다. 이를 확인해 보기 위해서 몇 가지 상황을 더 가정해보자.

Figure 16.	활동전위의 발생과 10ms일 때 1ms동안 15uA의 전기적 자극에 의한 막전위, gate 변수 그래프

① 25ms일 때 1ms동안 15[eqn: mu A]의 전류를 인가

이번에는 동일하게 1ms동안 15[eqn: mu A]의 전류를 인가하여 시뮬레이션하되 25ms 부근에서 전류를 인가하여 진행해보았다. 그 결과 Figure 17과 같은 그래프가 나타났다. 이 경우에는 처음(t=0일 때)에 막전위를 –70mV에서 –55mV로 증가시켰기 때문에 활동전위가 먼저 발생하고 그 이후 25ms 시점에서 전류를 가하여 막전위를 –68mV에서 -53mV로 증가시켰고 그에 따라 약 27mV까지 활동전위를 발생시켰다. gate variable m,h,n의 그래프는 전류가 인가된 부분에서 m과 n은 소폭 증가하고 h는 소폭 감소하는 모습을 보였고 이 외에는 앞서 관찰할 수 있었던 그래프가 동일하게 두 번 발생하는 것을 확인할 수 있었다.

이로써 동일한 크기의 전류를 동일한 시간동안 가했음에도 이렇게 다른 결과가 나타나는 것은 전류를 가했던 10ms를 지나는 이 시점이 불응기였기 때문이다. 불응기란 활동 전위 발생 직후 다음 활동전위가 발생하지 않는 수msec의 기간을 말한다[6]. 그렇기 때문에 이 기간이 지나고 25ms에 전류를 가했을 때는 활동전위는 발생하지만 10ms일 때는 활동전위가 발생하지 않는다.

Figure 17.	활동전위의 발생과 25ms일 때 1ms, 15uA의 전기적 자극에 의한 막전위, gate 변수 그래프

② 10ms일 때 2ms 동안 15[eqn: mu A]의 전류를 인가

먼저 앞선 시뮬레이션에서는 1ms동안 15[eqn: mu A]의 전류를 인가하였고 이번에는 2ms동안 전류의 크기는 동일하게 15[eqn: mu A]의 전류를 인가하여 시뮬레이션을 진행해보았다. 그 결과 Figure 18과 같은 그래프가 나타났다. 이 경우에는 처음(t=0일 때)에 막전위를 –70mV에서 –55mV로 증가시켰기 때문에 활동전위가 먼저 발생하고 그 이후 10ms 시점에서 전류를 가하여 2ms동안 막전위를 약 –58mV로 증가시켰고 그에 따라 약 23mV까지 활동전위를 발생시켰다. gate variable m,h,n의 그래프는 전류가 인가된 부분에서 m과 n은 소폭 증가하고 h는 소폭 감소하는 모습을 보였고 이 외에는 앞서 관찰할 수 있었던 그래프가 동일하게 두 번 발생하는 것을 확인할 수 있었다.

불응기에는 아무리 강한 자극을 주더라도 활동전위가 발생하지 않는 절대적 불응기와 상대적으로 강력한 자극을 주면 활동전위가 발생하는 상대적 불응기가 있다[5]. 이를 확인해보기 위에 전류 자극의 강도를 높여 2ms동안 자극이 지속되도록 한 것이다. 그 결과 활동전위는 발생하였기 때문이 10ms인 이 시점이 상대적으로 강력한 자극을 주면 활동전위가 발생하는 상대적 불응기임을 확인할 수 있다. 자극의 지속시간을 늘리는 방법을 통해서 활동전위를 발생시키는 더 강력한 자극을 만들어낼 수 있다는 것을 확인할 수 있다.

Figure 18.	활동전위의 발생과 10ms일 때 2ms, 15uA의 전기적 자극에 의한 막전위, gate 변수 그래프 (위) (아래) 표기 및 나누어서 설명

③ 10ms일 때 1ms 동안 30[eqn: mu A]의 전류를 인가

이번에는 1ms동안 전류의 크기를 더 크게 하여 30[eqn: mu A]의 전류를 인가하여 시뮬레이션을 진행해보았다. 그 결과 Figure 19과 같은 그래프가 나타났다. 이 경우에는 처음(t=0일 때)에 막전위를 –70mV에서 –55mV로 증가시켰기 때문에 활동전위가 먼저 발생하고 그 이후 10ms 시점에서 전류를 가하여 1ms동안 막전위를 약 –57mV로 증가시켰고 그에 따라 약 24mV까지 활동전위를 발생시켰다. gate variable m,h,n의 그래프는 전류가 인가된 부분에서 m과 n은 소폭 증가하고 h는 소폭 감소하는 모습을 보였고 이 외에는 앞서 관찰할 수 있었던 그래프가 동일하게 두 번 발생하는 것을 확인할 수 있었다.

앞선 시뮬레이션②를 통해 10ms가 상대적 불응기임을 확인하였고 자극의 지속시간을 늘리는 방법으로 활동전위를 발생시키는 더 강력한 자극을 만들어낼 수 있다는 것을 확인하였다. 이번에는 전류의 세기를 더 크게하여 시뮬레이션을 진행하였고 그 결과 시뮬레이션②와 동일하게 활동전위가 발생하는 것을 확인할 수 있다. 즉, 자극의 세기를 늘리는 방법을 통해서 활동전위를 발생시키는 더 강력한 자극을 만들어낼 수 있다는 것을 확인할 수 있었다.

Figure 19.	활동전위의 발생과 10ms일 때 1ms, 30uA의 전기적 자극에 의한 막전위, gate 변수 그래프 (위),... (아래)....혹은 (a) ... (b) ???

(Exercise 3) threshold 전압찾기

세포막에서 활동전위는 [eqn: Na ^{+}]통로의 positive feedback을 통해 모든 전압의존성  [eqn: Na ^{+}]통로를 활성화시킬만큼 초기의 막전위가 변화하지 않으면 발생하지 않는 실무율(all-or-nothing)을 따른다. 이때 필요한 막전위의 변화는 15mV에서 30mV의 급격한 변화이다. 일반적으로 threshold potenial은 –50mV~-55mV사이이며 다양한 요인에 의해 달라질 수 있다. 이를 알아보기 위해서 몇 가지 조건을 설정하여 시뮬레이션을 진행해보고자 한다.

① ip=0uA, vhold=-70mV, vstart=-70mV

이 경우에는 막전위를 voltage step increse시키지 않고 –70mV로 일정하게 유지할 때를 가정한다. 이에 따라 막전위는 큰 변화없이 약 –70mV로 일정하게 유지되며 이에 따라 gate-variable m,h,n도 0.07, 0.6 ,0.3으로 변함없이 유지된다. 이는 Figure 20의 그래프를 통해 확인할 수 있다.

Figure 20.	vstart=-70mV인 경우의 막전위와 gate variable그래프

② ip=0uA, vhold=-70mV, vstart=-67mV

이 경우에는 vstart를 –67mV로 하였기 때문에 처음 막전위가 –67mV이고 그 이후에 –70mV로 회복된다. 이는 활동전위를 일으키기에 충분하지 않았고 gate variable도 약간의 값의 변화는 있었지만 큰 변화없이 일정하게 유지되는 모습을 보였다.

Figure 21.	vstart=-67mV 경우의 막전위와 gate variable그래프

③ ip=0uA, vhold=-70mV, vstart=-64mV

이 경우에는 vstart를 –64mV로 하였기 때문에 처음 막전위가 –64mV이고 그 이후에 –70mV로 회복된다. 이 역시도 활동전위를 일으키기에 충분하지 않았고 gate variable도 약간의 값의 변화는 있었지만 큰 변화없이 일정하게 유지되는 모습을 보였다.

Figure 22.	vstart=-64mV일 경우의 막전위와 gate variable그래프

④ ip=0uA, vhold=-70mV, vstart=-63mV

이 경우 vstart를 –63mV로 하였기 때문에 처음 막전위가 –63mV이고 그 이후에 막전위는 급격하게 증가하여 3.8ms동안 약 +26mV에 도달하게 된다. gate variable에서도 [eqn: Na ^{+}]통로와 관련있는 m은 4.3ms 만에 0에서 1에 가까운 값까지 증가하게 되고, 반대로 h는 5ms동안 0.07로 매우 작은 값으로 감소하게 된다. 반면 [eqn: K ^{+}]통로와 관련있는 n은 5ms만에 약0.77로 증가하게 되며 서서히 감소하는 모습을 보인다. (Figure 23)

Figure 23.	vstart=-63mV일 경우의 막전위와 gate variable그래프

초기막전위 –70mV에서 7mV정도 증가시켰을 때 m의 값이 1에 가까운 값으로 증가한 것으로 보아 모든 전압의존성 [eqn: Na ^{+}]통로를 활성화시켰고 그에 따라 [eqn: Na ^{+}]는 세포안으로 유입되어 막전위가 +26mV까지 증가하는 탈분극이 일어날 수 있도록 한 것을 의미한다. 즉 본 프로젝트에서 모델링한 신경세포는 –70mV에서 –63mV로 빠르게 변화할 때 활동전위가 발생하며 이때 –63mV를 자극의 threshold라고 할 수 있다.

이러한 threshold 전위가 발생하는 원인에 대해서 고려해보자. Threshold 이하의 전압을 가해줄 경우를 가정하면 세포안으로 [eqn: Na ^{+}]가 이동하면서 흐르게 되는 전류가 증가하고 세포밖으로 [eqn: K ^{+}]와 다른 leakage ion들이 유출되면서 세포밖으로 흐르는 전류도 증가하게 된다. 그에 따라 발생하는 전류는 세포막을 불안정성하게 만드는데 이를 해소하기 위해 필요한 전위가 바로 threshold potential인 것이다[7].

(Exercise 4) Hyperpolarization

지금까지는 안정막 전위보다 높은 전압이 가해졌을 때 발생하는 탈분극과 그에 따라 이어지는 재분극에 의해 원래 상태로 회복되는 것을 보았다면 이번에는 그 반대 현상인 hyperpolarizing pulse에 따른 막전위와 gate variable들의 값들을 시뮬레이션하고자 한다.

vstart의 값을 –73mV에서 –88mV까지 –3mV씩 증분하여 이에 따른 membrane potential과 gate variable, 그리고 sodium, potassium conductance를 확인해보았다. 그 결과 –73mV,-76mV,-79mV,-82mV,-85mV,-88mV에서 모두 막전위와 gate variable은 초반부분의 미세한 변동이외에는 일정한 값을 유지하는 그래프가 나타났다. 즉, 이 모든 경우에서 action potential은 발생하지 않았다. 한편 sodium, potassium conductance의 경우 –88mV에 가까워 질수록 각각의 conductance는 미세하게 그 peak의 값이 커지는 것을 확인할 수 있었다.

	vstart=-73mV
	vstart=-76mV
	vstart=-79mV
	vstart=-82mV
	vstart=-85mV
	vstart=-88mV
Figure 24.	vstart의 값의 변화(-73mV~-88mV)에 따른 막전위, gate variable, ion conductance의 그래프

최종적으로 vstart=-90mV에서는 활동전위가 발생하였고 이에 따른 그래프는 Figure 25이다. 특이한 것은 impulse가 발생한 시점은 0ms인데 10ms 부근에서 활동전위가 발생하였으며 gate variable도 이 부근에서 값의 변동이 생긴다. 먼저 sodium channel과 관련한 m은 10ms가 되는 순간 급격하게 증가하여 급격하게 감소하는 모습을 보여준다. 한편 sodium channel의 inactivation과 관련된 h는 10ms부근에서부터 값이 감소하기 시작하여 sodium 통로의 문이 닫힘을 보여준다. 한편 potassium 통로와 관련된 n은 10ms 부근에서 값이 증가하여 완만하게 감소하는 형태를 보여주며 potassium을 세포막 밖으로 유출하여 막전위를 다시 본래 전위인 –70mV로 회복하는 것을 보여준다.

이러한 결과는 Figure 26에서 확인할 수 있듯이 sodium과 potassium 모두 10ms 부근에서 conductance의 변화가 발생하는 것을 보여주며 potassium conductance의 경우 sodium conductance에 비해 완만한 감소를 나타내는 것을 확인할 수 있다.

Figure 25.	vstart=-90mV일 경우의 막전위와 gate variable그래프

Figure 26.	vstart=-90mV일 경우의 Na+와 K+의 conductance

강한 세기의 hyperpolarizing stimulus도 자극이 끝나는 지점(trailing edge)에서 활동전위를 일으킬 수 있다. Hodgkin-Huxley는 이와 같은 현상을 anode break excitation 또는 rebound excitation라고 하였다. Hyperpolarizing pulse는 [eqn: K ^{+}]전도도를 감소시키고 [eqn: Na ^{+}] inactivation을 제거한다. [eqn: K ^{+}]전도도를 감소되면 hyperpolarizing current가 감소하고, [eqn: Na ^{+}] inactivation을 제거하게 되면 depolarizing current가 증가한다. [eqn: K ^{+}]통로의 움직임이 [eqn: Na ^{+}]통로의 움직임보다 느리기 때문에 일시적인 탈분극 현상이 발생하게 된다[8]. 이 시뮬레이션에서도 hyperpolarizing pulse가 시작되는 t=0이 아닌 10ms부근에서 활동전위가 발생하게 되는데 이 시점이 hyperpolarizing pulse가 끝나는 시점이기 때문에 위와 같은 그래프가 나타났다고 할 수 있다.

이러한 hyperpolarizing pulse는 depolarizing pulse와는 반대로 세포 내부의 electrode를 통해 전류가 전달될 때 조직 흥분의 활성화에 이용된다. 추가로 hyperpolarizing current는 rebound excitation mechanism을 통해 조직의 활성화로 이어질 수 있다[8].

CONCLUSIONS

본 프로젝트에서는 이러한 신경세포의 막의 구성, 물질운반, 전위, 흥분성 등의 특징으로 바탕으로 전기회로적으로 모델링하여 세포막에서 시간의 흐름에 따라 변화하는 막전위와 [eqn: Na ^{+} ,`K ^{+}]등의 conductance의 변화를 시뮬레이션해 보았다. 첫 번째로는 시뮬레이션 상에서 설계한 voltage step increase와 전류를 추가적으로 가함으로써 발생하는 변화를 관찰한 후 자극의 시점을 변화시키면서 10ms가 지나는 시점이 불응기라는 것을 확인할 수 있었다. 또한 인가한 전류의 세기와 지속시간을 2배로 늘리는 것을 통해 자극의 강도를 높였을 때 활동전위가 발생하는 것을 통해 이 시점이 상대적불응기라는 것과 자극의 강도를 높이는 방법에 대해서 알 수 있었다. 두 번째로는 시뮬레이션 상에서 설정한 전압 impulse를 안정상태 –70mV에서 3mV씩 증가시켜가며 활동전위가 어떤 조건 상에서 발생하는지 확인하였다. 그 결과 막전압이 –70mV에서 –63mV로 급격히 증가할 때 활동전위가 발생하는 것을 확인하였고 이를 모델링한 신경세포막의 역치 전위라고 할 수 있었다. 세 번째로는 두 번째와는 다르게 전압 impulse를 –70mV에서 반대방향으로 –3mV씩 증분해가며 활동전위가 발생하는 조건을 확인해보았다. 그 결과 막전위가 –70mV에서 –90mV로 급격히 증가하는 step impulse를 주었을 때 해당 impulse가 끝나는 시점에서 활동전위가 발생하는 rebound excitement가 나타나는 것을 확인할 수 있었다. 본 프로젝트에서 Hodgkin-Huxley model을 기반으로 신경세포막의 전위를 모델링하였는데 이 모델이 생리학적으로 설명하기는 어려운 모델이지만 실제 실험 데이터를 바탕으로 fitting을 진행하였고 시뮬레이션 결과가 실제 세포막과 유사한 결과를 나타냈기 때문에 simulation을 통해 얻어진 전위와 conductance의 값을 통해 세포막 사이의 ion의 이동을 통해 발생하는 활동전위와 여러 가지 전기자극을 통한 세포막 전위의 변화에 대해서 이해할 수 있었다.

REFERENCES

[1] John E.Hall, 의학생리학, 범문에듀케이션, 2012, p11~p26

[2] John E.Hall, 의학생리학, 범문에듀케이션, 2012, p63~77

[3] ‘Action potential and ion conductance’,memorangapp, viewed 2020.05.28

https://www.memorangapp.com/flashcards/94205/PHYS1+Resting+Membrane+Potential%2C+Action+%26+Graded+Potentials/

[4] ‘Electrical circuit model of the cell membrane’, researchgate, viewed 2020.05.28.

https://www.researchgate.net/figure/Electrical-circuit-model-of-the-cell-membrane-C-m-is-the-membrane-capacitance-R-m-is_fig9_327558169

[5] James Keener, Mathematical Physiology: I: Cellular Physiology, Springer Science & Business Media, 2010. 6. 4, p203

[6]‘불응기’, 생명과학대사전, 2020.05.31

https://terms.naver.com/entry.nhn?docId=427883&cid=60261&categoryId=60261

[7] ‘threshold potential’, wikipedia, 2020.05.31

https://en.wikipedia.org/wiki/Threshold_potential

[8] David Prutchi, John Wiley & Sons, Design and Development of Medical Electronic Instrumentation, 2005.1.28., p307

APPENDICES

(Exercise1)

HH.m

%filename HH.m

%numerical solution of the space-clamped Hodgkin-Huxley equations

clear all

clf

global check;

global t1p t2p ip; %parameters for the function izero(t)

in_HH

in_mhnv

for klok=1:klokmax

t=klok*dt;                      %note time

m=snew(m,alpham(v),betam(v),dt); %update m

h=snew(h,alphah(v),betah(v),dt); %update h

n=snew(n,alphan(v),betan(v),dt); %update n

gNa=gNabar*(m^3)*h;    %sodium conductance

gK =gKbar*(n^4);    %potassium conductance

g=gNa+gK+gLbar;         %total conductance

gE=gNa*ENa+gK*EK+gLbar*EL;         %gE=g*E

%save old value of v for checking purposes:

v_old=v;

%update v:

v=(v+(dt/C)*(gE+izero(t)))/(1+(dt/C)*g);

if(check)

E=gE/g;

chv=C*(v-v_old)/dt+g*(v-E)-izero(t)

end

%store results for future plotting:

mhn_plot(:,klok)=[m h n]';

v_plot(klok)=v;

t_plot(klok)=t;

end

subplot(2,1,1),plot(t_plot,v_plot)

title(subplot(2,1,1), 'Membrane Potential')

subplot(2,1,1), xlabel('t [ms]','Fontweight','bold'),ylabel('Membrane Potential[mV]','Fontweight','bold')

subplot(2,1,2),plot(t_plot,mhn_plot)

title(subplot(2,1,2), 'm,h,n의 그래프')

subplot(2,1,2), xlabel('t [ms]','Fontweight','bold'),ylabel('gate variable','Fontweight','bold')

in_HH.m

%filename: in_HH.m

%Units:

%voltage is in millivolts (mV)

%current is in microamperes (muA)

%length is in centimeters (cm)

%time is in milliseconds (ms)

%note that microfarads =

%microamperes*milliseconds/millivolt

%= muA*ms/mV

%

%initialize membrane parameters:

%membrane capacitance per unit area:

C=1.0      %(muF/cm^2)

%max possible Na+ conductance per unit area:

gNabar=120 %((muA/mV)/cm^2)

%max possible K+ conductance per unit area:

gKbar=36   %((muA/mV)/cm^2)

%leakage conductance per unit area:

gLbar=0.3  %((muA/mV)/cm^2)

%Na+ equilibrium potential:

ENa = 45   %(mV)

%K+ equilibrium potential:

EK = -82   %(mV)

%leakage channel reversal potential:

EL = -59   %(mV)

%

%initialize time step and experiment duration:

dt=0.1     %time step duration (ms)

tmax=50    %duration of experiment (ms)

%total number of time steps in the experiment:

klokmax=ceil(tmax/dt)

%

%initialize arrays that hold data for plotting:

mhn_plot=zeros(3,klokmax);

v_plot=zeros(1,klokmax);

t_plot=zeros(1,klokmax);

%

%initialize parameters that define the experiment:

%the neuron is at rest (v= -70 mV) prior to t=0;

%at t=0 a current shock is applied after which v= -55 mV;

%then a subsequent 15 muA current pulse of 1 ms duration

%is applied beginning at t=10 ms.

%voltage prior to t=0:

vhold=  -70  %(mV)

%voltage just after t=0:

vstart= -55%(mV)

%(change in v is result of current shock applied at t=0)

%

%initialize parameters of subsequent current pulse:

t1p=10       %starting time (ms)

t2p=11       %stopping time (ms)

ip=0        %current applied (muA)

%

%initialize checking parameter

check=1      %set check=1 to enable self-checking

%set check=0 to disable self-checking

in_mhnv.m

%filename: in_mhnv.m

%this script simulates holding the

%voltage constant for a long time

%prior to t=0 at the value v=vhold,

%and then applying a current shock

%at t=0 to step the voltage

%to a new level, v=vstart.

%The gating variables m,h,n

%cannot change suddenly and therefore

%remain constant during the current shock.

%

%simulate time prior to t=0:

v=vhold

%set m,h,n equal to their steady values

%under constant-v conditions:

m=alpham(v)/(alpham(v)+betam(v))

h=alphah(v)/(alphah(v)+betah(v))

n=alphan(v)/(alphan(v)+betan(v))

%

%now let voltage jump to its value

%just after t=0, without making

%any further change in m,h,n:

v=vstart

izero.m

function i=izero(t)

%filename: izero.m

global t1p t2p ip;

if((t1p<t)&(t<t2p))

i=ip;

else

i=0;

end

snew.m

function s=snew(s_old,alpha,beta,dt)

%filename: snew.m

global check;

s=(s_old+dt*alpha)/(1+dt*(alpha+beta));

if(check)

chsnew=(s-s_old)/dt-(alpha*(1-s)-beta*s)

end

alpham.m

function a=alpham(v)

%filename: alpham.m

theta=(v+45)/10;

if(theta==0)   %check for case that gives 0/0

a=1.0;  %in that case use L'Hospital's rule

else

a=1.0*theta/(1-exp(-theta));

end

alphah.m

function a=alphah(v)

%filename: alphah.m

theta = (v+70)/20;

a=0.07*exp(-theta);

alphan.m

function a=alphan(v)

%filename: alphan.m

theta=(v+60)/10;

if(theta==0)   %check for case that gives 0/0

a=0.1;  %in that case use L'Hospital's rule

else

a=0.1*theta/(1-exp(-theta));

end

betam.m

function b=betam(v)

%filename: betam.m

theta = (v+70)/18;

b=4.0*exp(-theta);

betah.m

function b=betah(v)

%filename: betah.m

theta = (v+40)/10;

b=1.0/(1+exp(-theta));

betan.m

function b=betan(v)

%filename: betan.m

theta = (v+70)/80;

b=0.125*exp(-theta);

(Exercise2)

in_HH.m _ t2p-t1p=2ms일 때

%filename: in_HH.m

%Units:

%voltage is in millivolts (mV)

%current is in microamperes (muA)

%length is in centimeters (cm)

%time is in milliseconds (ms)

%note that microfarads =

%microamperes*milliseconds/millivolt

%= muA*ms/mV

%

%initialize membrane parameters:

%membrane capacitance per unit area:

C=1.0      %(muF/cm^2)

%max possible Na+ conductance per unit area:

gNabar=120 %((muA/mV)/cm^2)

%max possible K+ conductance per unit area:

gKbar=36   %((muA/mV)/cm^2)

%leakage conductance per unit area:

gLbar=0.3  %((muA/mV)/cm^2)

%Na+ equilibrium potential:

ENa = 45   %(mV)

%K+ equilibrium potential:

EK = -82   %(mV)

%leakage channel reversal potential:

EL = -59   %(mV)

%

%initialize time step and experiment duration:

dt=0.1     %time step duration (ms)

tmax=50    %duration of experiment (ms)

%total number of time steps in the experiment:

klokmax=ceil(tmax/dt)

%

%initialize arrays that hold data for plotting:

mhn_plot=zeros(3,klokmax);

v_plot=zeros(1,klokmax);

t_plot=zeros(1,klokmax);

%

%initialize parameters that define the experiment:

%the neuron is at rest (v= -70 mV) prior to t=0;

%at t=0 a current shock is applied after which v= -55 mV;

%then a subsequent 15 muA current pulse of 1 ms duration

%is applied beginning at t=10 ms.

%voltage prior to t=0:

vhold=  -70  %(mV)

%voltage just after t=0:

vstart= -55  %(mV)

%(change in v is result of current shock applied at t=0)

%

%initialize parameters of subsequent current pulse:

t1p=10       %starting time (ms)

t2p=12       %stopping time (ms)

ip=15        %current applied (muA)

%

%initialize checking parameter

check=1      %set check=1 to enable self-checking

%set check=0 to disable self-checking

in_HH.m _ ip=30uA 일 때

%filename: in_HH.m

%Units:

%voltage is in millivolts (mV)

%current is in microamperes (muA)

%length is in centimeters (cm)

%time is in milliseconds (ms)

%note that microfarads =

%microamperes*milliseconds/millivolt

%= muA*ms/mV

%

%initialize membrane parameters:

%membrane capacitance per unit area:

C=1.0      %(muF/cm^2)

%max possible Na+ conductance per unit area:

gNabar=120 %((muA/mV)/cm^2)

%max possible K+ conductance per unit area:

gKbar=36   %((muA/mV)/cm^2)

%leakage conductance per unit area:

gLbar=0.3  %((muA/mV)/cm^2)

%Na+ equilibrium potential:

ENa = 45   %(mV)

%K+ equilibrium potential:

EK = -82   %(mV)

%leakage channel reversal potential:

EL = -59   %(mV)

%

%initialize time step and experiment duration:

dt=0.1     %time step duration (ms)

tmax=50    %duration of experiment (ms)

%total number of time steps in the experiment:

klokmax=ceil(tmax/dt)

%

%initialize arrays that hold data for plotting:

mhn_plot=zeros(3,klokmax);

v_plot=zeros(1,klokmax);

t_plot=zeros(1,klokmax);

%

%initialize parameters that define the experiment:

%the neuron is at rest (v= -70 mV) prior to t=0;

%at t=0 a current shock is applied after which v= -55 mV;

%then a subsequent 15 muA current pulse of 1 ms duration

%is applied beginning at t=10 ms.

%voltage prior to t=0:

vhold=  -70  %(mV)

%voltage just after t=0:

vstart= -55  %(mV)

%(change in v is result of current shock applied at t=0)

%

%initialize parameters of subsequent current pulse:

t1p=10       %starting time (ms)

t2p=11       %stopping time (ms)

ip=30        %current applied (muA)

%

%initialize checking parameter

check=1      %set check=1 to enable self-checking

%set check=0 to disable self-checking

(Exercise3)

in_HH.m

%filename: in_HH.m

%Units:

%voltage is in millivolts (mV)

%current is in microamperes (muA)

%length is in centimeters (cm)

%time is in milliseconds (ms)

%note that microfarads =

%microamperes*milliseconds/millivolt

%= muA*ms/mV

%

%initialize membrane parameters:

%membrane capacitance per unit area:

C=1.0      %(muF/cm^2)

%max possible Na+ conductance per unit area:

gNabar=120 %((muA/mV)/cm^2)

%max possible K+ conductance per unit area:

gKbar=36   %((muA/mV)/cm^2)

%leakage conductance per unit area:

gLbar=0.3  %((muA/mV)/cm^2)

%Na+ equilibrium potential:

ENa = 45   %(mV)

%K+ equilibrium potential:

EK = -82   %(mV)

%leakage channel reversal potential:

EL = -59   %(mV)

%

%initialize time step and experiment duration:

dt=0.1     %time step duration (ms)

tmax=50    %duration of experiment (ms)

%total number of time steps in the experiment:

klokmax=ceil(tmax/dt)

%

%initialize arrays that hold data for plotting:

mhn_plot=zeros(3,klokmax);

v_plot=zeros(1,klokmax);

t_plot=zeros(1,klokmax);

%

%initialize parameters that define the experiment:

%the neuron is at rest (v= -70 mV) prior to t=0;

%at t=0 a current shock is applied after which v= -55 mV;

%then a subsequent 15 muA current pulse of 1 ms duration

%is applied beginning at t=10 ms.

%voltage prior to t=0:

vhold=  -70  %(mV)

%voltage just after t=0:

vstart= -70  %(mV)

%(change in v is result of current shock applied at t=0)

%

%initialize parameters of subsequent current pulse:

t1p=10       %starting time (ms)

t2p=11       %stopping time (ms)

ip=0        %current applied (muA)

%

%initialize checking parameter

check=1      %set check=1 to enable self-checking

%set check=0 to disable self-checking

(Exercise4)

HH.m - gNa,gK그래프 그리기

%filename HH.m

%numerical solution of the space-clamped Hodgkin-Huxley equations

clear all

clf

global check;

global t1p t2p ip; %parameters for the function izero(t)

in_HH

in_mhnv

for klok=1:klokmax

t=klok*dt;                      %note time

m=snew(m,alpham(v),betam(v),dt); %update m

h=snew(h,alphah(v),betah(v),dt); %update h

n=snew(n,alphan(v),betan(v),dt); %update n

gNa=gNabar*(m^3)*h;    %sodium conductance

gK =gKbar*(n^4);    %potassium conductance

g=gNa+gK+gLbar;         %total conductance

gE=gNa*ENa+gK*EK+gLbar*EL;         %gE=g*E

%save old value of v for checking purposes:

v_old=v;

%update v:

v=(v+(dt/C)*(gE+izero(t)))/(1+(dt/C)*g);

if(check)

E=gE/g;

chv=C*(v-v_old)/dt+g*(v-E)-izero(t)

end

%store results for future plotting:

mhn_plot(:,klok)=[m h n]';

v_plot(klok)=v;

t_plot(klok)=t;

gNa_plot(klok)=gNa;

gK_plot(klok)=gK;

end

subplot(2,1,1)plot(t_plot, gNa_plot)

title(subplot(2,1,1), 'gNa')

subplot(2,1,1), xlabel('t [ms]','Fontweight','bold'),ylabel('gNa[(uA/mV)/cm^2]','Fontweight','bold')

subplot(2,1,2)plot(t_plot, gK_plot)

title(subplot(2,1,2), 'gK')

subplot(2,1,2), xlabel('t [ms]','Fontweight','bold'),ylabel('gK[(uA/mV)/cm^2]','Fontweight','bold')