Skip to content

Instantly share code, notes, and snippets.

@DaisukeMiyamoto

DaisukeMiyamoto/hodgkin-huxley.c

Last active June 10, 2017 02:08
Show Gist options
  • Save DaisukeMiyamoto/795e4d1122f845825054518706a4f007 to your computer and use it in GitHub Desktop.
Save DaisukeMiyamoto/795e4d1122f845825054518706a4f007 to your computer and use it in GitHub Desktop.
Download ZIP
calculate hodgkin huxley neuron model by python and C code based on NEURON simulator's method
Raw
hodgkin-huxley.c
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "mpi.h"
#ifdef KCOMPUTER
#include "fj_tool/fapp.h"
#endif
#ifdef USE_FLOAT
#define EXP(x) expf((x))
typedef float FLOAT;
#else
#define EXP(x) exp((x))
typedef double FLOAT;
#endif
//#define EXP(x) hoc_Exp((x))
double hoc_Exp(double x);
#define N_COMPARTMENT (1024 * 8 * 8 * 2)
static const FLOAT DT = 0.025; // [msec]
static const FLOAT CELSIUS = 6.3; // [degC]
static FLOAT v_trap (FLOAT x, FLOAT y) { return( (fabs(x/y) > 1e-6)? ( x/(EXP(x/y) - 1.0) ) : ( y*(1. - x/y/2.) ) ); }
static FLOAT calc_alpha_n (FLOAT v) { return( 0.01 * v_trap(-(v+55.), 10.) ); }
static FLOAT calc_beta_n (FLOAT v) { return( 0.125 * EXP( -(v+65.) / 80.) ); }
static FLOAT calc_alpha_m (FLOAT v) { return( 0.1 * v_trap(-(v+40.), 10.)); }
static FLOAT calc_beta_m (FLOAT v) { return( 4. * EXP( -(v+65) / 18.) ); }
static FLOAT calc_alpha_h (FLOAT v) { return( 0.07 * EXP( -(v+65) / 20.) ); }
static FLOAT calc_beta_h (FLOAT v) { return( 1. / (EXP( -(v+35) / 10.) + 1.) ); }
static FLOAT hh_v[N_COMPARTMENT]; // [mV]
static FLOAT hh_dv[N_COMPARTMENT];
static FLOAT hh_n[N_COMPARTMENT];
static FLOAT hh_m[N_COMPARTMENT];
static FLOAT hh_h[N_COMPARTMENT];
static FLOAT hh_cm[N_COMPARTMENT];
static FLOAT hh_cm_inv[N_COMPARTMENT];
static FLOAT hh_gk_max[N_COMPARTMENT];
static FLOAT hh_gna_max[N_COMPARTMENT];
static FLOAT hh_gm[N_COMPARTMENT];
const FLOAT hh_e_k = -77.0; // [mV]
const FLOAT hh_e_na = 50.0; // [mV]
const FLOAT hh_v_rest = -54.3; // [mV]
double hoc_Exp(double x){
if (x < -700.) {
return 0.;
}else if (x > 700) {
return exp(700.);
}
return exp(x);
}
static void initialize()
{
int i;
for(i=0; i<N_COMPARTMENT; i++)
{
hh_v[i] = -65.0; // [mV]
//hh_n[i] = 0.32;
//hh_m[i] = 0.05;
//hh_h[i] = 0.60;
hh_n[i] = calc_alpha_n(hh_v[i]) / (calc_alpha_n(hh_v[i]) + calc_beta_n(hh_v[i]));
hh_m[i] = calc_alpha_m(hh_v[i]) / (calc_alpha_m(hh_v[i]) + calc_beta_m(hh_v[i]));
hh_h[i] = calc_alpha_h(hh_v[i]) / (calc_alpha_h(hh_v[i]) + calc_beta_h(hh_v[i]));
hh_cm[i] = 1.0; // [muF/cm^2]
hh_cm_inv[i] = 1.0 / hh_cm[i]; // [cm^2/muF]
hh_gk_max[i] = 36.; // [mS/cm^2]
hh_gna_max[i] = 120.; // [mS/cm^2]
hh_gm[i] = 0.3; // [mS/cm^3]
}
return;
}
#define TABLE_SIZE 201
#define TABLE_MAX_V 100.0f
#define TABLE_MIN_V -100.0f
#ifdef TABLE_TYPE
FLOAT hh_table[TABLE_SIZE][6];
#define TABLE_N_TAU(x) hh_table[(x)][0]
#define TABLE_N_INF(x) hh_table[(x)][1]
#define TABLE_M_TAU(x) hh_table[(x)][2]
#define TABLE_M_INF(x) hh_table[(x)][3]
#define TABLE_H_TAU(x) hh_table[(x)][4]
#define TABLE_H_INF(x) hh_table[(x)][5]
#else
FLOAT hh_table[6][TABLE_SIZE];
#define TABLE_N_TAU(x) hh_table[0][(x)]
#define TABLE_N_INF(x) hh_table[1][(x)]
#define TABLE_M_TAU(x) hh_table[2][(x)]
#define TABLE_M_INF(x) hh_table[3][(x)]
#define TABLE_H_TAU(x) hh_table[4][(x)]
#define TABLE_H_INF(x) hh_table[5][(x)]
#endif
static void makeTable()
{
int i;
for(i=0; i<TABLE_SIZE; i++)
{
FLOAT v;
FLOAT a_n, a_m, a_h, b_n, b_m, b_h;
v = (TABLE_MAX_V - TABLE_MIN_V)/(FLOAT)(TABLE_SIZE-1) * i + TABLE_MIN_V;
a_n = calc_alpha_n(v);
a_m = calc_alpha_m(v);
a_h = calc_alpha_h(v);
b_n = calc_beta_n(v);
b_m = calc_beta_m(v);
b_h = calc_beta_h(v);
TABLE_N_TAU(i) = 1. / (a_n + b_n);
TABLE_N_INF(i) = a_n * TABLE_N_TAU(i);
TABLE_M_TAU(i) = 1. / (a_m + b_m);
TABLE_M_INF(i) = a_m * TABLE_M_TAU(i);
TABLE_H_TAU(i) = 1. / (a_h + b_h);
TABLE_H_INF(i) = a_h * TABLE_H_TAU(i);
//printf("%d : v(%.4f) n(%.4f %.4f) m(%.4f %.4f) h(%.4f %.4f)\n", i, v, TABLE_N_TAU(i), TABLE_N_INF(i), TABLE_M_TAU(i), TABLE_M_INF(i), TABLE_H_TAU(i), TABLE_H_INF(i) );
}
return;
}
static FLOAT calc_i_inj(FLOAT i_inj)
{
const FLOAT i_dc = 10.;
const FLOAT i_rand = 0.;
const FLOAT i_tau = 1.;
FLOAT i_inj_raw;
i_inj_raw = i_dc + ((float) rand() / RAND_MAX - 0.5) * 2. * i_rand;
i_inj += DT / i_tau * (i_inj_raw - i_inj);
return(i_inj);
}
typedef struct _memb{
FLOAT m;
FLOAT h;
FLOAT n;
FLOAT gk_max;
FLOAT gna_max;
FLOAT gm;
FLOAT v;
FLOAT cm_inv;
} Memb;
Memb memb[N_COMPARTMENT];
int hh_with_table(FLOAT stoptime)
{
unsigned int i, j;
unsigned int i_stop;
unsigned int v_i_array[N_COMPARTMENT];
FLOAT theta_array[N_COMPARTMENT];
const int inj_start = 50./DT;
const int inj_stop = 175./DT;
for(i=0,i_stop=stoptime/DT; i<i_stop; i++)
{
FLOAT i_inj;
if(i > inj_start && i < inj_stop){
i_inj = 10.0;
}else{
i_inj = 0.0;
}
//printf("%f %f %f %f\n", i*DT, i_inj, hh_v[0], hh_v[N_COMPARTMENT-1]);
#pragma omp parallel
{
#ifdef KCOMPUTER
fapp_start("prepare", 4, 2);
#endif
#pragma omp for
for(j=0; j<N_COMPARTMENT; j++)
{
v_i_array[j] = (int)(hh_v[j] - TABLE_MIN_V);
theta_array[j] = (hh_v[j] - TABLE_MIN_V) - (FLOAT)v_i_array[j];
}
#pragma omp for
for(j=0; j<N_COMPARTMENT; j++)
{
if(!(v_i_array[j] >= TABLE_SIZE || v_i_array[j]<0) ){
;
}else if(v_i_array[j] >= TABLE_SIZE){
v_i_array[j]=TABLE_SIZE-1; theta_array[j]=1.0;
}else if(v_i_array[j] < 0){
v_i_array[j]=0; theta_array[j]=0.0;
}
}
#ifdef KCOMPUTER
fapp_stop("prepare", 4, 2);
#endif
#ifdef KCOMPUTER
fapp_start("state", 2, 2);
#endif
#pragma omp for
for(j=0; j<N_COMPARTMENT; j++)
{
FLOAT tau_n, n_inf, tau_m, m_inf, tau_h, h_inf;
unsigned int v_i = v_i_array[j];
FLOAT theta = theta_array[j];
tau_n = TABLE_N_TAU(v_i);
n_inf = TABLE_N_INF(v_i);
tau_m = TABLE_M_TAU(v_i);
m_inf = TABLE_M_INF(v_i);
tau_h = TABLE_H_TAU(v_i);
h_inf = TABLE_H_INF(v_i);
tau_n = (tau_n + theta * (TABLE_N_TAU(v_i+1) - tau_n));
tau_m = (tau_m + theta * (TABLE_M_TAU(v_i+1) - tau_m));
tau_h = (tau_h + theta * (TABLE_H_TAU(v_i+1) - tau_h));
n_inf = n_inf + theta * (TABLE_N_INF(v_i+1) - n_inf) - hh_n[j];
m_inf = m_inf + theta * (TABLE_M_INF(v_i+1) - m_inf) - hh_m[j];
h_inf = h_inf + theta * (TABLE_H_INF(v_i+1) - h_inf) - hh_h[j];
hh_n[j] += (1.0f - EXP(-DT/tau_n)) * n_inf;
hh_m[j] += (1.0f - EXP(-DT/tau_m)) * m_inf;
hh_h[j] += (1.0f - EXP(-DT/tau_h)) * h_inf;
hh_v[j] += DT * hh_cm_inv[0] * (hh_gk_max[0] * hh_n[j] * hh_n[j] * hh_n[j] * hh_n[j] * (hh_e_k - hh_v[j])
+ hh_gna_max[0] * hh_m[j] * hh_m[j] * hh_m[j] * hh_h[j] * (hh_e_na - hh_v[j])
+ hh_gm[0] * (hh_v_rest - hh_v[j]));
}
#ifdef KCOMPUTER
fapp_stop("state", 2, 2);
#endif
}
}
return(0);
}
int main(int argc, char **argv)
{
//printf("Hodgkin-Huxley equation\n");
int myid;
FLOAT stoptime = 100;
MPI_Init(&argc, &argv);
MPI_Comm_rank(MPI_COMM_WORLD, &myid);
initialize();
makeTable();
printf (" [Hodgkin-Huxley Conditions]\n");
printf (" * dt = %f, nstep = %d\n", DT, (int)(stoptime/(FLOAT)DT));
printf (" * Total Compartments = %d\n", N_COMPARTMENT);
printf (" * Type : %s\n", (sizeof(FLOAT)==sizeof(double))?"double":"float");
#ifdef KCOMPUTER
fapp_start("calc", 1, 1);
#endif
printf("start (%d)\n", myid);
hh_with_table(stoptime);
printf("finished (%d)\n", myid);
#ifdef KCOMPUTER
fapp_stop("calc", 1, 1);
#endif
MPI_Finalize();
return(0);
}
Raw
hodgkin-huxley.py
# -*- coding: utf-8 -*-
"""
Created on Thu Sep 17 01:29:38 2015
@author: nebula
"""
import math
import matplotlib.pyplot as plt
class HodgkinHuxley():
LINESTYLE = ['-', ':', '--', '-.']
I_INJ = 5 # [nA]
DT = 0.025 # [msec]
def __init__(self):
self.table = {}
self.table_max_v = 100.0
self.table_min_v = -100.0
self.table_size = 201
self.i_inj = []
self.cm = 1.0 # [muF/cm^2]
self.e_k = -77.0 # [mV]
self.e_na = 50.0 # [mV]
self.gk_max = 36.0 # [mS/cm^2]
self.gna_max = 120.0 # [mS/cm^2]
self.gm = 0.3 # [mS/cm^3]
self.v_rest = -54.3 # [mV]
self.v = -65.0 # [mV]
self.n = self._alpha_n(self.v) / (self._alpha_n(self.v) + self._beta_n(self.v))
self.m = self._alpha_m(self.v) / (self._alpha_m(self.v) + self._beta_m(self.v))
self.h = self._alpha_h(self.v) / (self._alpha_h(self.v) + self._beta_h(self.v))
self.makeTable()
def _vtrap(self, x, y):
x = float(x)
y = float(y)
if math.fabs(x / y) > 1e-6:
return x / (math.exp(x/y) - 1.0)
else:
return y * (1.0 - x/y/2.)
def _alpha_n(self, v):
return 0.01 * self._vtrap(-(v+55.), 10.)
def _beta_n(self, v):
return 0.125 * math.exp(-(v+65.) / 80.)
def _alpha_m(self, v):
return 0.1 * self._vtrap(-(v+40.), 10.)
def _beta_m(self, v):
return 4. * math.exp(-(v+65.) / 18.)
def _alpha_h(self, v):
return 0.07 * math.exp(-(v+65.) / 20.)
def _beta_h(self, v):
return 1.0 / (math.exp(-(v+35.)/10.) + 1.0)
def _i2v(self, i):
return (self.table_max_v - self.table_min_v) / (self.table_size -1) * float(i) + self.table_min_v
def makeTable(self):
self.table['n_tau'] = []
self.table['n_inf'] = []
self.table['m_tau'] = []
self.table['m_inf'] = []
self.table['h_tau'] = []
self.table['h_inf'] = []
for i in range(self.table_size):
v = self._i2v(i)
a_n = self._alpha_n(v)
b_n = self._beta_n(v)
a_m = self._alpha_m(v)
b_m = self._beta_m(v)
a_h = self._alpha_h(v)
b_h = self._beta_h(v)
self.table['n_tau'].append(1.0 / (a_n + b_n))
self.table['n_inf'].append(a_n / (a_n + b_n))
self.table['m_tau'].append(1.0 / (a_m + b_m))
self.table['m_inf'].append(a_m / (a_m + b_m))
self.table['h_tau'].append(1.0 / (a_h + b_h))
self.table['h_inf'].append(a_h / (a_h + b_h))
def plotDict(self, dct):
i = 0
rng = [self._i2v(tmp_i) for tmp_i in range(self.table_size)]
for k,v in dct.items():
plt.plot(rng, v, linestyle=(self.LINESTYLE[i % len(self.LINESTYLE)]), linewidth=2, label=k)
i += 1
plt.xlabel('V [mV]')
plt.ylabel('Y')
plt.legend()
plt.show()
def showTable(self):
st = ' V :'
for k in self.table.keys():
for i in range(12 - len(k)):
st += '-'
st += k + '-'
st += '\n'
for i in range(0, self.table_size, 10):
st += '%10.4f :' % self._i2v(float(i))
for val in self.table.values():
st += ' %12.8f' % val[i]
st += '\n'
print st
def _i_inj(self):
t = self.t
if t > 5:
return self.I_INJ
else:
return 0.0
def _v2i(self, v):
i = int(v - self.table_min_v)
theta = (v - self.table_min_v) - float(i)
if(i >= self.table_size-2):
i = self.table_size-2
theta = 1.0
if(i < 0):
i = 0
theta = 0.0
#print i,theta
return i, theta
def calc_n(self, v, n):
table = self.table
i, theta = self._v2i(v)
n_tau = table['n_tau'][i] + theta * (table['n_tau'][i+1] - table['n_tau'][i])
n_inf = table['n_inf'][i] + theta * (table['n_inf'][i+1] - table['n_inf'][i])
n += (1.0 - math.exp(-self.DT/n_tau)) * (n_inf - n)
self.n = n
return n
def calc_m(self, v, m):
table = self.table
i, theta = self._v2i(v)
m_tau = table['m_tau'][i] + theta * (table['m_tau'][i+1] - table['m_tau'][i])
m_inf = table['m_inf'][i] + theta * (table['m_inf'][i+1] - table['m_inf'][i])
m += (1.0 - math.exp(-self.DT/m_tau)) * (m_inf - m)
self.m = m
return m
def calc_h(self, v, h):
table = self.table
i, theta = self._v2i(v)
h_tau = table['h_tau'][i] + theta * (table['h_tau'][i+1] - table['h_tau'][i])
h_inf = table['h_inf'][i] + theta * (table['h_inf'][i+1] - table['h_inf'][i])
h += (1.0 - math.exp(-self.DT/h_tau)) * (h_inf - h)
self.h = h
return h
def calc_v(self, v, n, m, h):
v += self.DT / self.cm \
* (self.gk_max * (n**4) * (self.e_k - v) \
+ self.gna_max * (m**3) * h * (self.e_na - v) \
+ self.gm * (self.v_rest - v) + self._i_inj())
self.v = v
return v
if __name__ == '__main__':
hh = HodgkinHuxley()
hh.showTable()
hh.plotDict(hh.table)
record = {}
record['t'] = [0]
record['v'] = [hh.v]
record['n'] = [hh.n]
record['m'] = [hh.m]
record['h'] = [hh.h]
for i in range(1000):
record['t'].append(record['t'][-1] + hh.DT)
hh.t = record['t'][-1]
record['n'].append(hh.calc_n(record['v'][-1], record['n'][-1]))
record['m'].append(hh.calc_m(record['v'][-1], record['m'][-1]))
record['h'].append(hh.calc_h(record['v'][-1], record['h'][-1]))
record['v'].append(hh.calc_v(record['v'][-1], record['n'][-1], record['m'][-1], record['h'][-1]))
plt.plot(record['t'], record['v'])
plt.ylim(-80, 80)
plt.xlabel('t [msec]')
plt.ylabel('V [mV]')
plt.show()
plt.plot(record['t'], record['m'], label='m')
plt.plot(record['t'], record['h'], label='h')
plt.plot(record['t'], record['n'], label='n')
n4 = [n**4 for n in record['n']]
m3h = [m**3*h for m,h in zip(record['m'], record['h'])]
plt.plot(record['t'], m3h, label='m^3 * h')
plt.plot(record['t'], n4, label='n^4')
plt.legend(loc=1)
plt.xlabel('t [msec]')
plt.show()
#plt.legend(loc=1)
#plt.xlabel('t [msec]')
#plt.show()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment