2017-03-16 1 views
3

파이썬에서 두 개의 뉴런 네트워크를 시뮬레이트하려고합니다. 그것은 각 뉴런에 대해 별도의 방정식을 작성하는 것만 큼 간단하지만 코드를 좀 더 일반화하여 방정식을 반복해서 다시 작성하지 않고도 뉴런의 수를 늘릴 수 있기를 바랍니다.파이썬 : 행렬로 작은 뉴럴 네트워크 작성하기

enter image description here

기본적으로, I는 전압 방정식에서 마지막 항을 통해 함께 결합 된 두 지킨 - 헉슬리 뉴런을 가지고 두 개의 신경망 방정식은 다음이다. 그래서 제가하고 싶은 것은 그러한 방식으로 코드를 작성하여 네트워크를 쉽게 확장 할 수 있도록하는 것입니다. 이를 위해, 뉴런 전압 [V1, V2]에 대한 벡터 V를 생성하고 벡터 X를 생성합니다. 여기서 X는 게이팅 변수 m, h 및 n을 모델링합니다. 그래서 저는 X = [[m1, h1, n1], [m2, h2, n2]]를 가질 것입니다. 그러나 현재 코드는 스파이크를 생성하지 않고 오히려 전압이 무한대로 불어 난 것처럼 보입니다. 이것은 게이팅 변수 X에 문제가 있음을 나타냅니다. 게이팅 변수 m, h 및 n은 항상 0과 1 사이 여야하므로 게이팅 변수가 1에 도달하고 거기에 머물러있는 것처럼 보입니다 쪽으로. 나는 그들이 1에 머물러야하는 원인이 무엇인지 확신하지 못한다. 코드가 실행 중이고 오류가 발생하지 않는다. 내가 의도적으로 나중에 때문에 방정식에 확률 성을 추가하고 싶습니다 때문에 내 미분 방정식을 통합하는 odeint를 사용하고 있지 않다

import scipy as sp 
import numpy as np 
import pylab as plt 

NN=2 #Number of Neurons in Model 

dt=0.01 
T = sp.arange(0.0, 1000.0, dt) 
nt = len(T) # total number of time steps 

    # Constants 
gNa = 120.0 # maximum conducances, in mS/cm^2 
gK = 36.0 
gL = 0.3 
ENa = 50.0 # Nernst reversal potentials, in mV 
EK = -77 
EL = -54.387 

#Coupling Terms 
Vr = 20 
w = 1 
e11 = e22 = 0 
e12 = e21 = 0.1 
E = np.array([[e11, e12], [e21, e22]]) 

#Gating Variable Transition Rates 
def alpham(V): return (0.1*V+4.0)/(1.0 - sp.exp(-0.1*V-4.0)) 
def betam(V): return 4.0*sp.exp(-(V+65.0)/18.0) 
def alphah(V): return 0.07*sp.exp(-(V+65.0)/20.0) 
def betah(V): return 1.0/(1.0 + sp.exp(-0.1*V-3.5)) 
def alphan(V): return (0.01*V+0.55)/(1.0 - sp.exp(-0.1*V-5.5)) 
def betan(V): return 0.125*sp.exp(-(V+65.0)/80.0) 
def psp(V,s): return ((5*(1-s))/(1+sp.exp(-(V+3)/8)))-s 

#Current Terms 
def I_Na(V,x): return gNa * (x[:,0]**3) * x[:,1] * (V - ENa) #x0=m, x1=h, x2=n 
def I_K(V,x): return gK * (x[:,2]**4) * (V - EK) 
def I_L(V): return gL * (V - EL) 
def I_inj(t): return 10.0 

#Initial Conditions 
V = np.zeros((nt,NN)) #Voltage vector 
X = np.zeros((nt,NN,3)) #Gating Variables m,h,n (NN neurons x 3 gating variables) 
S = np.zeros((nt,NN)) #Coupling term 

dmdt = np.zeros((nt,NN)) 
dhdt = np.zeros((nt,NN)) 
dndt = np.zeros((nt,NN)) 

V[0,:] = -65.0 
X[0,:,0] = alpham(V[0,:])/(alpham(V[0,:])+betam(V[0,:])) #m 
X[0,:,1] = alphah(V[0,:])/(alphah(V[0,:])+betah(V[0,:])) #h 
X[0,:,2] = alphan(V[0,:])/(alphan(V[0,:])+betan(V[0,:])) #n 

alef = 5.0/(1+sp.exp(-(V[0,:]+3)/8.0)) 
S[0,:] = alef/(alef+1) 

dmdt[0,:] = alpham(V[0,:])*(1-X[0,:,0])-betam(V[0,:])*X[0,:,0] 
dhdt[0,:] = alphah(V[0,:])*(1-X[0,:,1])-betah(V[0,:])*X[0,:,1] 
dndt[0,:] = alphan(V[0,:])*(1-X[0,:,2])-betan(V[0,:])*X[0,:,2] 

#Euler-Maruyama Integration 
for i in xrange(1,nt): 
    V[i,:]= V[i-1,:]+dt*(I_inj(i-1)-I_Na(V[i-1,:],X[i-1,:])-I_K(V[i-1,:],X[i-1,:])-I_L(V[i-1,:]))+dt*((Vr-V[i-1,:])/w * np.dot(E,S[i-1,:]))  

    #Gating Variable 
    dmdt[i,:] = dmdt[i-1,:] + alpham(V[i-1,:])*(1-X[i-1,:,0])-betam(V[i-1,:])*X[i-1,:,0] 
    dhdt[i,:] = dhdt[i-1,:] + alphah(V[i-1,:])*(1-X[i-1,:,1])-betah(V[i-1,:])*X[i-1,:,1] 
    dndt[i,:] = dndt[i-1,:] + alphan(V[i-1,:])*(1-X[i-1,:,2])-betan(V[i-1,:])*X[i-1,:,2] 
    z = np.array([dmdt[i-1,:],dhdt[i-1,:],dndt[i-1,:]]).T 

    #Gating Variable Constraints (0<m,h,n<1) 
    X[i,:,0] = max(0,min(X[i,:,0].all(),1)) 
    X[i,:,1] = max(0,min(X[i,:,1].all(),1)) 
    X[i,:,2] = max(0,min(X[i,:,2].all(),1)) 

    #Update Gating Variables 
    X[i,:,:]= X[i-1,:,:]+dt*(z) 

    #Coupling Term 
    S[i,:] = S[i-1,:]+dt*psp(V[i-i,:],S[i-1,:]) 

V1 = V[:,0] 
V2 = V[:,1] 

plt.plot(T,V1, 'red') 
plt.plot(T,V2, 'blue') 

plt.show() 

위의 오일러 방법을 사용하고 싶습니다. 어쨌든 누군가가이 코드를 수정하여 예상되는 스파이크 동작이 발생하는 것을 파악할 수 있다면 환상적 일 것입니다. 고맙습니다!

+0

생물 물리 모델 인 HH 모델에 대해 말하면 '신경 네트워크'라고 표시하는 것을 망설이게됩니다. –

+0

@ juanpa.arrivillaga 좋습니다. 내가 바 꾸었습니다. – Brenton

+0

먼저 코드의 개별 조각을 먼저 작동시킨 다음 단계별로 복잡성을 추가하는 것이 좋습니다. 예를 들어 수동 전류가 하나의 뉴런에 작용하는 것을 고려한 다음 한 뉴런에 대해 스파이 킹을 얻은 다음 두 번째 뉴런의 수동 전류를 추가 한 다음 두 가지 모두에 대해 스파이크 동작을 추가합니다. 문제는 방정식 또는 통합 코드와 관련이있을 수 있으며 단계별로 작업하는 것이 디버깅을 쉽게 처리 할 수 ​​있어야합니다. – Justas

답변

0

당신이

\dot(x)_j = f(x_j) + \sum_j C_ij g(x_j, x_i)

을 가지고 있기 때문에 당신은 또한 같은 당신의 타는 편을 쓸 수 있습니다. x_j는 벡터입니다 (예 : (v, m, h, n). 나는이 개선 될 수있다 확신 n_oscillators

def f(x): 
    return np.array([x[1, :], -x[0, :]]) 


def g(xi, xj): 
    return xi[0] - xj[0] 

def rhs(x, t, c): 
    coupling = np.array([sum(cij*g(xi,xj) for cij, xj in zip(ci, x.T)) 
              for ci, xi in zip(c, x.T)]) 
    coupling = np.outer(np.arange(2), coupling) # coupling in x'' 

    return f(x) + coupling 

x0 = np.random.random(size=(2, 3)) 
>>> array([[ 0.74386362, 0.85799588, 0.70501992], 
    [ 0.65903405, 0.41575505, 0.93166973]]) 

def ring(n): 
    c = np.eye(n, k=1) + np.eye(n, k=-1) 
    c[0, -1] = 1 
    c[-1, 0] = 1 
    return c 

c = ring(3) 
x1 = rhs(x0, 0, c) 
>>> array([[ 0.65903405, 0.41575505, 0.93166973], 
    [-0.81915217, -0.59088767, -0.89683958]]) 

,

x.shape이 n_onsite_variables입니다 :

이 오실레이터 예

는 출발점이 될 수있다. 특히 좀 더 일반적인 커플 링을 원한다면.

1

입력 전류 및 컨덕턴스를 확인하십시오. 현대화 된 형식주의에 따라 코드를 작성하면 dm/dt = (m_inf - m)/tau라는 더 쉽게 사용할 수 있습니다. 특히, 게이팅 변수 통합은 작동하지 않습니다. 당신은 그들을 제대로 업데이트하지 않습니다. timestep 수학이 누락되었는지 확인하십시오.

+0

도움을 주셔서 감사합니다! dmdt, dhdt 및 dndt에 대해 사용한 방정식이 올바르게 업데이트되지 않았습니다. 이제 해결되었습니다! – Brenton

관련 문제