Below is a script with a single compartment with a K channel after Hodgkin and Huxley 1952 (squid demo).
When the simulation is run repeatedly to obtain K channel conductance at different voltage steps, the value of Gk is retained after reinit. In each run, this value left from the previous run gradually falls to baseline.
Only if the model is explicitly scheduled, as in squid demo, does the Gk value get reset properly.
# hh_script.py ---
#
# Filename: hh_script.py
# Description:
# Author: Subhasis Ray
# Created: Sat Feb 17 11:36:01 2024 (+0530)
# Last-Updated: Sat Feb 17 14:31:10 2024 (+0530)
# By: Subhasis Ray
#
# Code:
"""This script is for testing Hodgkin-Huxley type K+ channel"""
# Usual imports
import numpy as np
import matplotlib.pyplot as plt
import moose
# Functions to setup voltage clamp
def create_voltage_clamp(modelpath, datapath, compartment):
"""Creates a voltage clamp object under `modelpath` and a table
under `datapath` to record the command voltage.
Returns the `moose.PulseGen` that gives the command value to the
voltage clamp.
"""
path = f'{modelpath}/vclamp'
if moose.exists(path):
print(f'{path}: Object already exists')
vclamp = moose.VClamp(path)
command = moose.PulseGen(f'{modelpath}/command')
command_tab = moose.Table(f'{datapath}/command')
return vclamp, command, command_tab
vclamp = moose.VClamp(path)
# The voltage clamp's output is `currentOut` which will be
# injected into the compartment
moose.connect(vclamp, 'currentOut', compartment, 'injectMsg')
# The voltage clamp object senses the voltage of the compartment
moose.connect(compartment, 'VmOut', vclamp, 'sensedIn')
command = moose.PulseGen(f'{modelpath}/command')
# Connect the output of the command pulse to the command input of the voltage clamp circuit
moose.connect(command, 'output', vclamp, 'commandIn')
# Also setup a table to record the command voltage of the VClamp directly
command_tab = moose.Table(f'{datapath}/command')
moose.connect(command_tab, 'requestOut', command, 'getOutputValue')
# compartment.dt is the integration time step for the compartment
# `tau` is the time constant of the RC filter in the circuit.
# 5 times the integration timestep value is a good starting point for tau
vclamp.tau = 5 * compartment.dt
# `gain` is the proportional gain of the PID controller. `Cm/dt` is a good value
vclamp.gain = compartment.Cm / compartment.dt
# `ti` is the integral time of the PID controller, `dt` is a good value
vclamp.ti = compartment.dt
# `td` is the derivative time of the PID controller. We can set it to 0
return vclamp, command, command_tab
def set_command_timecourse(command, base, delay, level):
"""Set up an existing pulse generator `command` to output `base`
as initial value and `level` after `delay` time
"""
command.baseLevel = base
command.firstDelay = delay
command.secondDelay = 1e9 # Never stop
command.firstWidth = 1e9 # Never stop
command.firstLevel = level
# Create containers
data = moose.Neutral('/data')
model = moose.Neutral('/model')
elec = moose.Neutral('/electronics')
# Create a compartment representing the squid giant axon
axon = moose.Compartment(f'{model.path}/axon')
# Create a Hodgkin-Huxley-type channel and connect it to the compartment
kchan = moose.HHChannel(f'{axon.path}/K')
moose.connect(kchan, 'channel', axon, 'channel')
# moose.showmsg(kchan.path)
kchan.Xpower = 4
kchan.Ek = -12.0 # mV with respect to resting Vm
# Gbar is maximum conductance.
# The voltage dependent conductance value is computed and
# stored in the field Gk, which we record in a table below
kchan.Gbar = 36.0 # mS/cm^2 -
gate = moose.HHGate(f'{kchan.path}/gateX')
# gate.useInterpolation = True # use a lookup table for alpha and beta
alpha_params = [0.1, -0.01, -1.0, -10.0, -10.0]
beta_params = [0.125, 0, 0, 0, 80.0]
vdivs = 150
vmin = -30.0
vmax = 120.0
# Note that `+` operator with lists as operands concatenates them
params = alpha_params + beta_params + [vdivs, vmin, vmax]
gate.setupAlpha(params)
# Insert voltage clamp circuitry
vclamp, command, command_tab = create_voltage_clamp(elec.path, '/data', axon)
gK_tab = moose.Table(f'{data.path}/gK')
moose.connect(gK_tab, 'requestOut', kchan, 'getGk')
im_tab = moose.Table(f'{data.path}/Im')
moose.connect(im_tab, 'requestOut', axon, 'getIm')
axon.Em = 0 # Hodgkin and Huxley used resting voltage as 0
axon.initVm = 0
axon.Cm = 1.0 # uF/cm^2
axon.Rm = 1 / 0.3 # G_leak is 0.3 mS/cm^2
# The voltage steps used in HH paper, relative to resting Vm
commands = np.array([6, 10, 19, 26, 32, 38, 51, 63, 76, 88, 100, 109], dtype=float)
fig, axes = plt.subplots(nrows=3, ncols=1, sharex='all')
runtime = 50.0
# ############################################################
## DEBUG START
# simdt = 1e-5
# plotdt = 1e-2
# moose.setClock(0, simdt)
# moose.setClock(1, simdt)
# moose.setClock(2, simdt)
# moose.setClock(3, plotdt)
# moose.useClock(0, f'{model.path}/#[TYPE=Compartment]', 'init')
# moose.useClock(0, f'{elec.path}/##', 'process')
# moose.useClock(1, f'{model.path}/#[TYPE=Compartment]', 'process')
# moose.useClock(2, f'{axon.path}/#[TYPE=HHChannel]', 'process')
# moose.useClock(3, f'{data.path}/#[TYPE=Table]', 'process')
## DEBUG END
# ############################################################
for level in commands:
# baseline is set to equilibrium potential(0.0)
# delay is 10 ms, start the voltage change at this time
# command level is relative to resting
set_command_timecourse(command, 0.0, 10.0, level)
moose.reinit()
print(f'At start of {level} mV: {kchan.Gk} mS/cm^2')
moose.start(runtime)
print(f'At end of: {level} mV: {kchan.Gk} mS/cm^2')
gK = gK_tab.vector
t_gK = np.arange(len(gK)) * gK_tab.dt
axes[0].plot(t_gK, gK, label=f'{level} mV')
im = im_tab.vector
t_im = np.arange(len(im)) * im_tab.dt
axes[1].plot(t_im, im)
v_command = command_tab.vector
t_command = np.arange(len(v_command)) * command_tab.dt
axes[2].plot(t_command, v_command)
axes[0].set_ylabel('Conductivity (mS/cm^2)')
axes[1].set_ylabel('Im (uA/cm^2)')
axes[1].set_ylim(0, 5000.0)
axes[2].set_xlabel('Time (ms)')
axes[2].set_ylabel('Command voltage (mV)')
fig.set_size_inches(4, 6)
fig.savefig('hhchannel_scheduling_issue.png')
plt.show()
# import sys
# sys.exit(0)
#
# hh_script.py ends here
Below is a script with a single compartment with a K channel after Hodgkin and Huxley 1952 (squid demo). When the simulation is run repeatedly to obtain K channel conductance at different voltage steps, the value of Gk is retained after
reinit. In each run, this value left from the previous run gradually falls to baseline.Only if the model is explicitly scheduled, as in squid demo, does the Gk value get reset properly.