Simulating (semi-)morphological dynamics in the Hodgkin-Huxley model

Support
7

Ah, that explains it! Good to know, and yes I agree that it might avoid confusion if this were specified in the documentation. :+1:

Sure, I will post the full current version, with plots and all (maybe I made some errors there, too?)

Also, another issue I can’t seem to avoid is that my voltage doesn’t cross -60 mV in the dynamics. I wonder if I forgot to provide some input and whether this is related to my issue of no dendritic influence on the dynamics?

# Neuron parameters:

defaultclock.dt = 0.01 * ms
v_0 = -65 * mV

# Morphology and specifications:

morpho = Soma(40 * um)
morpho.axon = Cylinder(diameter = 1*um, length = 300*um, n = 100)
morpho.L = Cylinder(diameter = 1*um, length = 150*um, n = 50) # Left dendrite.
morpho.R = Cylinder(diameter = 1*um, length= 150*um, n = 50) # Right dendrite

# Conductances:

# Conductances:
# Removed the other fixed parameters in order for them to remain unbounded.
g_L_const = 0.3 * msiemens * meter**-2

# Reversal potentials:

E_L = -54 * mV
E_K = -77 * mV
E_Na = 50 * mV

#____________________________________________________

# Neuron models & input equation(s):

# HH:

Ho_Hux = Equations(
'''
Im = ((g_Na * m**3 * h * (E_Na - v)) + (g_K * n**4 * (E_K - v)) + (g_L * (E_L - v))) : amp / meter**2

dm/dt = alpha_m * (1. - m) - beta_m * m : 1

dn/dt = alpha_n * (1. - n) - beta_n * n : 1

dh/dt = alpha_h * (1. - h) - beta_h * h : 1

alpha_n = (0.01 * (v / mV + 55) / (1 - exp(-(v / mV + 55) / 10))) / ms : Hz
beta_n = (0.125 * exp(-(v / mV + 65) / 80)) / ms : Hz

alpha_m = (0.1 * (v / mV + 40) / (1 - exp(-(v / mV + 40) / 10))) / ms : Hz
beta_m = (4.0 * exp(-0.0556 * (v / mV + 65))) / ms : Hz

alpha_h = (0.07 * exp(-0.05 * (v / mV + 65))) / ms : Hz
beta_h = (1 / (1 + exp(-(0.1) * (v / mV + 35)))) / ms : Hz

I : amp (point current)
g_K : siemens * meter**-2
g_Na : siemens * meter**-2
g_L : siemens * meter**-2

'''
)


# Instantiate neuron(s):

MorphNeuron = SpatialNeuron(morphology=morpho, model=Ho_Hux, 
                       Cm=1*uF/cm**2, Ri=200*ohm*cm, threshold='v > 0*mV', refractory='v > -10*mV', reset= 'v = -60*mV', method='exponential_euler')
 
 # Initialize group parameters: 
                      
MorphNeuron.v = v_0
MorphNeuron.I = 0 * amp

# Synapses and stimulation (later):


# Monitor(s):

somaMon = StateMonitor(MorphNeuron,  ['v', 'm', 'n', 'h'], record=True)
mon_L = StateMonitor(MorphNeuron.L, ['v', 'm', 'n', 'h'], record=True)
mon_R = StateMonitor(MorphNeuron.R, ['v', 'm', 'n', 'h'], record=True)

#____________________________________________________
# Running simulation and plotting:

# Current injections:

run(10 * ms)

MorphNeuron.I[morpho.R[50*um]] = 0.02 * nA  # injecting in the right dendrite

run(5 * ms)

MorphNeuron.I = 0.002 * amp
MorphNeuron.I[morpho.L[100*um]] = 0.01 * nA 
#MorphNeuron.I[0] = 0.2 * amp

# No influence?
MorphNeuron.L[10*um:30*um].g_Na = 500 * g_L_const
MorphNeuron.L[30*um:150*um].g_Na = 600 * g_L_const
MorphNeuron.R[30*um:60*um].g_Na = 300 * g_L_const
#MorphNeuron.L[60*um:230*um].g_Na = 700 * g_L_const # out of bounds actually?

run(50 * ms, report='text')

# Set (bigger) plot size:
plt.rcParams['figure.figsize'] = [15, 11]

# Plots:

subplot(4, 1, 1)
plot(somaMon.t/ms, somaMon[0].v/mV, label='Soma', color='k')
plot(mon_L.t/ms, mon_L[morpho.L[50*um]].v/mV, linestyle='dotted', label='Left dendrite', color='r')
plot(mon_R.t/ms, mon_R[morpho.R[50*um]].v/mV, linestyle='dashed', label='Right dendrite', color='b')
xlabel('Time (in ms)')
ylabel('v (in mV)')
plt.legend()

subplot(4, 1, 2)
plot(somaMon.t/ms, somaMon.m[0], label='m (Na Channels), Soma', color='b')
plot(somaMon.t/ms, somaMon.n[0], label='n (K Channels), Soma', color='r')
plot(somaMon.t/ms, somaMon.h[0], label='h (K Channels), Soma', color='k')
xlabel('Time (in ms)')
ylabel('Opening probability for gating variables')
plt.legend()

subplot(4, 1, 3)
plot(mon_L.t/ms, mon_L[morpho.L[100 * um]].m, label='m (Na Channels), LD', color='b') 
plot(mon_L.t/ms, mon_L[morpho.L[100 * um]].n, label='n (K Channels), LD', color='r')
plot(mon_L.t/ms, mon_L[morpho.L[100 * um]].h, label='h (K Channels), LD', color='k')
xlabel('Time (in ms)')
ylabel('Opening probability for gating variables')
plt.legend()

subplot(4, 1, 4)
plot(mon_R.t/ms, mon_R[morpho.R[100 * um]].m, label='m (Na Channels), RD', color='b')
plot(mon_R.t/ms, mon_R[morpho.R[100 * um]].n, label='n (K Channels), RD', color='r')
plot(mon_R.t/ms, mon_R[morpho.R[100 * um]].h, label='h (K Channels), RD', color='k')
xlabel('Time (in ms)')
ylabel('Opening probability for gating variables')
plt.legend()

show()