Hodgkin Huxley PymoNNto Implementation
The following code creates a Neuron with multiple Hodgkin Huxley segments. Over time the input current of the first compartment changes, which creates different spike patterns flowing through the segments.
# /Examples_Classical/Hodgkin_Huxley.py
from PymoNNto import *
import scipy as sp
import matplotlib.pyplot as plt
#multi compartment version of
#https://hodgkin-huxley-tutorial.readthedocs.io/en/latest/_static/Hodgkin%20Huxley.html
#with a simplified version of the cable theory to connect compartments
class HH_main(Behavior):
def alpha_m(self, V):#Channel gating kinetics. Functions of membrane voltage
return 0.1*(V+40.0)/(1.0 - np.exp(-(V+40.0) / 10.0))
def beta_m(self, V): #Channel gating kinetics. Functions of membrane voltage
return 4.0*np.exp(-(V+65.0) / 18.0)
def alpha_h(self, V): #Channel gating kinetics. Functions of membrane voltage
return 0.07*np.exp(-(V+65.0) / 20.0)
def beta_h(self, V): #Channel gating kinetics. Functions of membrane voltage
return 1.0/(1.0 + np.exp(-(V+35.0) / 10.0))
def alpha_n(self, V): #Channel gating kinetics. Functions of membrane voltage
return 0.01*(V+55.0)/(1.0 - np.exp(-(V+55.0) / 10.0))
def beta_n(self, V): #Channel gating kinetics. Functions of membrane voltage
return 0.125*np.exp(-(V+65) / 80.0)
def I_Na(self, V, m, h): #Membrane current (in uA/cm^2) Sodium (Na = element name)
return self.g_Na * m**3 * h * (V - self.E_Na)
def I_K(self, V, n): #Membrane current (in uA/cm^2) Potassium (K = element name)
return self.g_K * n**4 * (V - self.E_K)
def I_L(self, V):#Membrane current (in uA/cm^2) Leak
return self.g_L * (V - self.E_L)
def I_inj(self, i): #External Current
return 5*(i*self.dt>100) - 5*(i*self.dt>200) + 10*(i*self.dt>300) - 10*(i*self.dt>400)
def initialize(self, n):
self.set_parameters_as_variables(self)
#C_m = 1.0 # membrane capacitance, in uF/cm^2
#g_Na = 120.0 # Sodium (Na) maximum conductances, in mS/cm^2
#g_K = 36.0 # Postassium (K) maximum conductances, in mS/cm^2
#g_L = 0.3 # Leak maximum conductances, in mS/cm^2
#E_Na = 50.0 # Sodium (Na) Nernst reversal potentials, in mV
#E_K = -77.0 # Postassium (K) Nernst reversal potentials, in mV
#E_L = -54.387 # Leak Nernst reversal potentials, in mV
#b_r = 0.13 #Resistance between blocks in ohms between blocks
self.dt = 0.01
n.v = np.array([n.vector()*self.v for i in range(self.blocks)])
n.h = np.array([n.vector()*self.h for i in range(self.blocks)])
n.m = np.array([n.vector()*self.m for i in range(self.blocks)])
n.n = np.array([n.vector()*self.n for i in range(self.blocks)])
n.I = np.array([n.vector() for i in range(self.blocks)])
def iteration(self, neurons):
v, h, m, n = neurons.v.copy(), neurons.h.copy(), neurons.m.copy(), neurons.n.copy()
I = neurons.I*0
#external current on first block of all neurons in group
I[0, :] = self.I_inj(neurons.iteration)
# Warning: this is a very simplified way to connect the compartments.
# The version is broadly based on the cable theory with a simple resistance between the compartments.
# It is NOT physically accurate and only there to demonstrate some form of spike propagation.
# Do not use this transition function for scientific research without further testing!
# current input to next block
for i in range(self.blocks-1):
I[i + 1] += ((v[i]-v[i+1])/self.b_r) * self.dt
# current input from previous block
for i in range(self.blocks-1):
I[i] += ((v[i + 1]-v[i])/self.b_r) * self.dt
##############################################################
neurons.v += ((I - self.I_Na(v, m, h) - self.I_K(v, n) - self.I_L(v)) / self.C_m) * self.dt
neurons.m += (self.alpha_m(v)*(1.0-m) - self.beta_m(v)*m) * self.dt
neurons.h += (self.alpha_h(v)*(1.0-h) - self.beta_h(v)*h) * self.dt
neurons.n += (self.alpha_n(v)*(1.0-n) - self.beta_n(v)*n) * self.dt
My_Network = Network()
N_e = NeuronGroup(net=My_Network, tag='excitatory_neurons', size=1, behavior={
1: HH_main(blocks=20, b_r=0.15, v=-65, m=0.05, h=0.6, n=0.32, C_m=1.0, g_Na=120.0, g_K=36.0, g_L=0.3, E_Na=50.0, E_K=-77.0, E_L=-54.387),
9: Recorder(['v','m','h','n'], tag='my_recorder')
})
SynapseGroup(net=My_Network, src=N_e, dst=N_e, tag='GLUTAMATE')
My_Network.initialize()
My_Network.simulate_iterations(50000, measure_block_time=True)
data = My_Network['v', 0, 'np'][:, :, 0]#numpy array of first recorder \ only the forst neuron and all compartments
for d in range(data.shape[1]):
plt.plot(data[:, d]-10*d)
plt.show()
plt.imshow(data.transpose(), cmap='gray', aspect='auto', interpolation='none')
plt.show()