|
| 1 | +import scipy as sp |
| 2 | +import numpy as np |
| 3 | +import pylab as plt |
| 4 | +from scipy.integrate import odeint |
| 5 | + |
| 6 | +class HodgkinHuxley(): |
| 7 | + """Full Hodgkin-Huxley Model implemented in Python""" |
| 8 | + |
| 9 | + C_m = 1.0 |
| 10 | + """membrane capacitance, in uF/cm^2""" |
| 11 | + |
| 12 | + g_Na = 120.0 |
| 13 | + """Sodium (Na) maximum conductances, in mS/cm^2""" |
| 14 | + |
| 15 | + g_K = 36.0 |
| 16 | + """Postassium (K) maximum conductances, in mS/cm^2""" |
| 17 | + |
| 18 | + g_L = 0.3 |
| 19 | + """Leak maximum conductances, in mS/cm^2""" |
| 20 | + |
| 21 | + E_Na = 50.0 |
| 22 | + """Sodium (Na) Nernst reversal potentials, in mV""" |
| 23 | + |
| 24 | + E_K = -77.0 |
| 25 | + """Postassium (K) Nernst reversal potentials, in mV""" |
| 26 | + |
| 27 | + E_L = -54.387 |
| 28 | + """Leak Nernst reversal potentials, in mV""" |
| 29 | + |
| 30 | + t = np.arange(0.0, 450.0, 0.01) |
| 31 | + """ The time to integrate over """ |
| 32 | + |
| 33 | + def alpha_m(self, V): |
| 34 | + """Channel gating kinetics. Functions of membrane voltage""" |
| 35 | + return 0.1*(V+40.0)/(1.0 - np.exp(-(V+40.0) / 10.0)) |
| 36 | + |
| 37 | + def beta_m(self, V): |
| 38 | + """Channel gating kinetics. Functions of membrane voltage""" |
| 39 | + return 4.0*np.exp(-(V+65.0) / 18.0) |
| 40 | + |
| 41 | + def alpha_h(self, V): |
| 42 | + """Channel gating kinetics. Functions of membrane voltage""" |
| 43 | + return 0.07*np.exp(-(V+65.0) / 20.0) |
| 44 | + |
| 45 | + def beta_h(self, V): |
| 46 | + """Channel gating kinetics. Functions of membrane voltage""" |
| 47 | + return 1.0/(1.0 + np.exp(-(V+35.0) / 10.0)) |
| 48 | + |
| 49 | + def alpha_n(self, V): |
| 50 | + """Channel gating kinetics. Functions of membrane voltage""" |
| 51 | + return 0.01*(V+55.0)/(1.0 - np.exp(-(V+55.0) / 10.0)) |
| 52 | + |
| 53 | + def beta_n(self, V): |
| 54 | + """Channel gating kinetics. Functions of membrane voltage""" |
| 55 | + return 0.125*np.exp(-(V+65) / 80.0) |
| 56 | + |
| 57 | + def I_Na(self, V, m, h): |
| 58 | + """ |
| 59 | + Membrane current (in uA/cm^2) |
| 60 | + Sodium (Na = element name) |
| 61 | +
|
| 62 | + | :param V: |
| 63 | + | :param m: |
| 64 | + | :param h: |
| 65 | + | :return: |
| 66 | + """ |
| 67 | + return self.g_Na * m**3 * h * (V - self.E_Na) |
| 68 | + |
| 69 | + def I_K(self, V, n): |
| 70 | + """ |
| 71 | + Membrane current (in uA/cm^2) |
| 72 | + Potassium (K = element name) |
| 73 | +
|
| 74 | + | :param V: |
| 75 | + | :param h: |
| 76 | + | :return: |
| 77 | + """ |
| 78 | + return self.g_K * n**4 * (V - self.E_K) |
| 79 | + # Leak |
| 80 | + def I_L(self, V): |
| 81 | + """ |
| 82 | + Membrane current (in uA/cm^2) |
| 83 | + Leak |
| 84 | +
|
| 85 | + | :param V: |
| 86 | + | :param h: |
| 87 | + | :return: |
| 88 | + """ |
| 89 | + return self.g_L * (V - self.E_L) |
| 90 | + |
| 91 | + def I_inj(self, t): |
| 92 | + """ |
| 93 | + External Current |
| 94 | +
|
| 95 | + | :param t: time |
| 96 | + | :return: step up to 10 uA/cm^2 at t>100 |
| 97 | + | step down to 0 uA/cm^2 at t>200 |
| 98 | + | step up to 35 uA/cm^2 at t>300 |
| 99 | + | step down to 0 uA/cm^2 at t>400 |
| 100 | + """ |
| 101 | + return 10*(t>100) - 10*(t>200) + 35*(t>300) - 35*(t>400) |
| 102 | + |
| 103 | + @staticmethod |
| 104 | + def dALLdt(X, t, self): |
| 105 | + """ |
| 106 | + Integrate |
| 107 | +
|
| 108 | + | :param X: |
| 109 | + | :param t: |
| 110 | + | :return: calculate membrane potential & activation variables |
| 111 | + """ |
| 112 | + V, m, h, n = X |
| 113 | + |
| 114 | + dVdt = (self.I_inj(t) - self.I_Na(V, m, h) - self.I_K(V, n) - self.I_L(V)) / self.C_m |
| 115 | + dmdt = self.alpha_m(V)*(1.0-m) - self.beta_m(V)*m |
| 116 | + dhdt = self.alpha_h(V)*(1.0-h) - self.beta_h(V)*h |
| 117 | + dndt = self.alpha_n(V)*(1.0-n) - self.beta_n(V)*n |
| 118 | + return dVdt, dmdt, dhdt, dndt |
| 119 | + |
| 120 | + def Main(self): |
| 121 | + """ |
| 122 | + Main demo for the Hodgkin Huxley neuron model |
| 123 | + """ |
| 124 | + |
| 125 | + X = odeint(self.dALLdt, [-65, 0.05, 0.6, 0.32], self.t, args=(self,)) |
| 126 | + V = X[:,0] |
| 127 | + m = X[:,1] |
| 128 | + h = X[:,2] |
| 129 | + n = X[:,3] |
| 130 | + ina = self.I_Na(V, m, h) |
| 131 | + ik = self.I_K(V, n) |
| 132 | + il = self.I_L(V) |
| 133 | + |
| 134 | + plt.figure() |
| 135 | + |
| 136 | + ax1 = plt.subplot(4,1,1) |
| 137 | + plt.title('Hodgkin-Huxley Neuron') |
| 138 | + plt.plot(self.t, V, 'k') |
| 139 | + plt.ylabel('V (mV)') |
| 140 | + |
| 141 | + plt.subplot(4,1,2, sharex = ax1) |
| 142 | + plt.plot(self.t, ina, 'c', label='$I_{Na}$') |
| 143 | + plt.plot(self.t, ik, 'y', label='$I_{K}$') |
| 144 | + plt.plot(self.t, il, 'm', label='$I_{L}$') |
| 145 | + plt.ylabel('Current') |
| 146 | + plt.legend() |
| 147 | + |
| 148 | + plt.subplot(4,1,3, sharex = ax1) |
| 149 | + plt.plot(self.t, m, 'r', label='m') |
| 150 | + plt.plot(self.t, h, 'g', label='h') |
| 151 | + plt.plot(self.t, n, 'b', label='n') |
| 152 | + plt.ylabel('Gating Value') |
| 153 | + plt.legend() |
| 154 | + |
| 155 | + plt.subplot(4,1,4, sharex = ax1) |
| 156 | + i_inj_values = [self.I_inj(t) for t in self.t] |
| 157 | + plt.plot(self.t, i_inj_values, 'k') |
| 158 | + plt.xlabel('t (ms)') |
| 159 | + plt.ylabel('$I_{inj}$ ($\\mu{A}/cm^2$)') |
| 160 | + plt.ylim(-1, 40) |
| 161 | + |
| 162 | + plt.tight_layout() |
| 163 | + plt.show() |
| 164 | + |
| 165 | +if __name__ == '__main__': |
| 166 | + runner = HodgkinHuxley() |
| 167 | + runner.Main() |
0 commit comments