Amor Natura

 이 논문은 허드킨-헉슬리 모델을 비롯한 여러 neuron spiking model들과 마찬가지로 neuron의 spiking 패턴을 미분방정식으로 만들어 모델링한 내용을 담고 있습니다. 그러나 내용이 짧아 왜 그렇게 되는지는 모르겠습니다. 물론 Izhikevich가 만든 모델을 이용해 다른 여러 현상들을 보는게 저의 목적이기 때문에 깊이 이해하고 있지는 않습니다. 그러나 간단히 그가 만든 모델이 어떤 것이고 무엇을 설명할 수 있는지 정리할 필요가 있습니다.

1. The Model

---

이 논문에 소개된 모델은 다른 여러 모델들보다 훨씬 간결하고 계산이 쉽습니다. Hodgkin-Huxley Model의 경우 연립 미분방정식이 4개로 이루어져 있습니다만 이 모델은 2개의 연립 미분방정식으로 뉴런의 발화(firing)패턴을 거의 비슷하게 만들 수 있습니다. 

                (1)

                                              (2)

그리고 부가적으로 뉴런이 발화한 후에 변수를 바꾸어 줍니다. 

여기서 v와 u는 차원이 없는 변수들이며 a, b, c, d는 차원이 없는 parameter 입니다. 

v와 u에 미분형태는 시간에 대한 미분(d/dt)입니다. v는 막전위(membrane potential)을 의미하며 u는 막회복변수(membrane recovery variable)을 의미합니다. 논문에서는 변수 u는 칼륨 이온의 화성화와 나트륨 이온의 비활성화에 대해 설명하고 있으며 u는 v를 억제하는 피드백을 준다고 하는군요. 

(1)의 I는 시냅스에서의 전류나 뉴런에 주어진 직류 전류를 의미합니다.

이렇게 변수들의 설명은 끝났고 이제 parameter에 대해서 정리를 해보도록 하죠.

a는 u에 대한 time scale을 의미한다고 합니다. 작은 값일 수록 u가 회복되는 시간이 오래 걸리는 것 같죠? (2)식에서 보면 u'과 u는 a만큼의 관계가 있으니 말입니다. 

실제로 뉴런의 모델링에선 0.02로 둡니다.

b는 논문에서 u가 얼마나 민감한지를 막전위 v의 subthreshold fluctuation으로 설명하는 parameter라고 합니다. b값이 커질수록 u와 v의 관계가 커지며 이 결과로 subthreshold 진동과 low-threshold spiking dynamics가 가능하다고 합니다. 보통 뉴런의 경우 0.2 로 둡니다.

c는 막전위인 v가 발화(firing, spiking)후 되돌아가는 초기값을 의미합니다. c는 -65로 둡니다. 이해를 돕기 위해 단위를 mV로 생각할 수 있습니다.

d는 회복변위인 u가 발화 후 되돌아가는 초기값을 의미합니다. 보통 2로 둡니다. 

2. Different Types of Neurons

---

  이제 두 개의 연립미분방정식으로 만들어진 Izhikevich의 모델을 이용해 여러 뉴런들을 계산해보기로 합니다. 우선 뉴런의 종류부터 알아야 하는데요. 뉴런은 크게 흥분성(Excitatory Cortical Cell)과 억제성(Inhibitory Cortical Cell)로 나누어 집니다. 흥분성 뉴런의 경우 자신이 받은 전류들을 모아 자신과 연결된 뉴런들에게 전달해주는 역할을 합니다. 그러나 억제성 뉴런의 경우에는 자신이 받은 전류들을 자신과 연결된 뉴런들에게 전달해주지 않습니다.

 이제 이 두가지 뉴런의 spiking type을 Izhikevich모델을 통해 구현해보겠습니다.

 ① RS(regular spiking) Neurons

 RS 뉴런은 대뇌피질(cortex)에서 가장 많이 발견되는 전형적인 뉴런입니다. 파라미터는 다음과 같습니다.

a = 0.02,   b = 0.2,  c = -65,   d = 8

이를 C언어를 이용해 미분방정식을 풀어보면 다음과 같은 결과를 얻습니다.

x축은 시간축으로 단위는 ms 이며 y축은 전압(mV)을 의미합니다. 미분방정식을 푼 method는 뒤에 소개하도록 하겠습니다.

② CH(chttering) Neurons

CH 뉴런은 spike가 일어나는 부근에서 아주 가까운 시간에 burst가 일어나는 패턴을 나타냅니다. 각 burst들 사이의 frequency는 40Hz정도 된다고 합니다.

파라미터는 다음과 같습니다.

a = 0.02,   b = 0.2,   c = -50,   d = 2

마찬가지로 C언어를 이용해 미분방정식을 풀어보면 실제와 비슷한 결과를 얻습니다.

①FS(Fast Spiking) Neurons

FS 뉴런은 RS와 달리 순응과정(adaptation, slowing down)과정이 없이 계속 빠르게 Spiking하는 뉴런입니다. 

파라미터는 다음과 같습니다.

a = 0.1,   b = 0.2,   c = -65,   d = 2

마찬가지로 미분방정식을 컴퓨터를 통해 풀어보면 다음과 같은 결과를 얻습니다.

이외에도 논문에는 몇 가지 다른 흥분성, 억제성 뉴런이 소개되어 있으며 이외에도 thalamo-cortical neuron 도 소개되어 있습니다. 그러나 주어진 parameter를 이용해 계산을 해보니 실제와 비슷한 결과가 나오지 않았습니다... 이유를 모르겠네요... 아무튼 흥분성 뉴런과 억제성 뉴런의 대표인 RS와 FS가 잘 나온것으로도 충분히 좋은 모델이라고 생각됩니다.

다른 모델에 비해 빠르고 간단하기 때문에 많은 수의 뉴런을 시뮬레이션하기에 적합한 것 같습니다.

(출처 :: http://www.izhikevich.org/publications/spikes.htm  )

3. Pulse-Coupled Implementation

Excitatory neuron은 RS Neuron으로, Inhibitory Neuron은 FS Neuron으로 가정합니다. 또한 Ex/In Neuron의 비율은 80:20으로 정합니다. 이 비율은 Thalamus 에서 Ex/In의 비율이 80:20 임에서 착안한 것 같습니다. 

각각의 RS, FS neuron의 parameter는 다음과 같이 정해집니다. 논문에 이렇게 시뮬레이션했다고 나오는군요..

흥분성 뉴런(RS, CH neuron)	억제성 뉴런 (FS Neuron)

위의 표에서  는  사이의 uniformly random variable 입니다. 흥분성 뉴런의 경우  값이 0이면 RS Neuron, 1이면 CH Neuron입니다. 

이외에도 주목할만한 변수와 미분방정식을 푼 method를 소개해야 하는데요. 논문에서는 이 Fully Connected network 에서 각 Neuron들이 연결되어 있는 Weight를 다음과 같이 주었습니다. 위에서와 마찬가지로 random number를 이용해 만듭니다.

 이외에도 각 뉴런에 주어지는 전류 역시 다음과 같이 주어집니다. 전류는 지금까지의 parameter와 달리 Gaussian Random Distribution을 이용합니다. Uniform하지 않죠.

 또 미분방정식을 푸는 해법이 좀 독특합니다. Euler method를 이용해 미분방정식을 푸는데요. 이상하게도 v에 대해선 dt를 줄여 2번 게산합니다. u에 대해서도 dt를 줄여서 계산하거나 Euler Method 대신 Runge-Kutta 2nd Methdod를 쓰면 원래 패턴이 잘 재현되지 않네요. 이에 대한 논문도 교수님께서 보여주신 적이 있었습니다. 그러나 이 방법대로 하면 여러 패턴들이 잘 나오니 상당히 유용하다는 것은 분명합니다.

1. 
     v=v+0.5*(0.04*v.^2+5*v+140-u+I); % step 0.5 ms
     v=v+0.5*(0.04*v.^2+5*v+140-u+I); % for numerical
     u=u+a.*(b.*v-u);                 % stability
   
   

위의 코드는 http://www.izhikevich.org/publications/net.m 에서 보실 수 있으며 MATLAB으로 만들어졌습니다. C로 코딩하면서 이와 같은 방법으로 미분방정식을 풀었으며 논문과 비슷한 결과를 얻을 수 있었습니다.

위의 그림은 제가 C언어로 다시 만들어서 얻은 결과이고, 아래 그림은 논문에 있는 결과입니다.  모두 raster plot 이며 세로축은 각 뉴런의 number, 가로축은 시간축을 나타냅니다. 발화(firing)할 때마다 점을 찍으면 위의 그림처럼 나타나게 되지요. 논문의 저자는 이 raster plot 에서 처음 보이는 진한 라인이 알파파(10Hz)이고 뒤에 감마파(40Hz)가 만들어지는 것을 볼 수 있다면서 포유류의 뇌파도 이 모델을 통해 만들 수 있다고 이야기를 합니다. 

자세한 내용은 그의 홈페이지(http://www.izhikevich.org )에서 확인하실 수 있습니다.

이 글은 스프링노트에서 작성되었습니다.

표랑 일부 서식이 깨져서 나오네요.. 스프링노트에서 더 깔끔하게 보실 수 있습니다.

수리물리에서 자주 등장하는 Symbol 중 가장 기초적인 것은 Kroneker-Delta(크로네커-델타)와 Levi-Civita Symbol(레비 시비타 심볼)이 있습니다.

우선 약속부터 보면,

(1) Kroneker-Delta

크로네커 델타의 경우 인덱스의 두개가 같으면 1이고 다르면 0이 됩니다.

(2) Levi-Civita

레비-시비타의 경우 permutation이라는 단어가 나옵니다. 순열이라는 뜻인데요. 이는 1,2,3 이 순서대로 나오는 것을 의미합니다.

따라서 123 혹은 231, 312를 의미합니다. 숫자의 순서가 바뀌지 않으면서 배열되어야 하지요.

이와 반대의 경우는 permutation이 아닌 경우(반대의 경우) 입니다. 즉 213, 213, 321 인 경우이지요.

3개의 숫자를 세 번 써서 만들 수 있는 경우의 수 중에 이 6가지를 제외한 다른 방법은 모두 같은 숫자가 2개 이상 들어가게 됩니다.

따라서 이런 경우에는 값이 0이 됩니다.

이런 약속은 왜 쓰냐구요? 크로네커 델타의 경우는 다양하게 쓰입니다. 예를 들자면 직교성을 살펴볼 때 사용할 수 있습니다.

레비-시비타는 벡터의 외적(corss product)의 계산을 빠르게 하기 위해 쓰입니다. 외적은 다음과같이 계산될 수 있습니다.

이 식을 천천히 뜯어보도록 하겠습니다. 우선 i=1 일 때는 x축 성분값을 구할 수 있습니다. 그럼 j와 k는 2, 3 또는 3,2 가 되겠지요. 숫자가 서로 2,2 또는 3,3일 경우 입실론은 0이 됩니다. 그럼 계산은 다음처럼 됩니다.

마찬가지로 계산을 따라가면 다음과 같은 계산을 얻을 수 있습니다. 이는 외적의 계산과 똑같지요.

이외에 응용으로는 다음과 같은 방법이 있으며 계산이 복잡한 문제를 푸는데 큰 도움을 줍니다.

이상입니다.

이 글은 스프링노트에서 작성되었습니다.