Description of problem
I am trying to implement an HH model.
I have compiled some channel models to build it.
I have verified that they are correct, that they produce the expected traces in VC.
Whenever I add a channel model with a gating particle (p) equation of the pattern:
d p/dt = (pinf-p/tau)
The resulting model will fail to equilibrate in CC.
gating particle values assume large negative values although they should be in a range of 0 to1 always.
Minimal code to reproduce problem
# %% code
from brian2 import *
# data frame support and `tidy` plotting
import pandas as pd
import numpy as np
from scipy import signal
import seaborn as sns
# Assign some channels
paramNaT = dict()
# m particle
paramNaT["Aam"] = 0.182/mvolt # A paramenter of m alpha rate function
paramNaT["Abm"] = 0.124/mvolt # A parameter of betam rate function
paramNaT["Vam"] = -35*mvolt # V1/2 of alpha m
paramNaT["Vbm"] = -35*mvolt # V1/2 of beta m
paramNaT["kam"] = 9*mvolt # slope of alpha m
paramNaT["kbm"] = 9*mvolt # slope of beta m
# h particle
paramNaT["Vhinf"] = -65*mvolt # hinf Boltzmann params
paramNaT["khinf"] = 6.2*mvolt #
paramNaT["Aah"] = 0.024 / mvolt # A paramenter of h alpha rate function
paramNaT["Abh"] = 0.0091/mvolt # A parameter of h beta rate function
paramNaT["Vah"] = -50*mvolt # V1/2 of alpha h
paramNaT["Vbh"] = -75*mvolt # V1/2 of beta h
paramNaT["kah"] = 5*mvolt # slope of alpha h
paramNaT["kbh"] = 5*mvolt # slope of beta h
eqINaT = '''
#v : volt # comment out if needed for CC experis
# Current equation
INaT = gNaT * m**3*h*(v-ENa) : amp
# activation particle equations
dm/dt = alpham*(1-m)-betam*m : 1
alpham = Aam*kam/(exprel(-(v-Vam)/kam))/ms : Hz
# old wrong exprel version
#betam = (-Abm*(v-Vbm))/(1-exp((v-Vbm)/kbm))/ms : Hz
# correct
betam = Abm*kbm/exprel((v-Vbm)/kbm)/ms : Hz
# inactivation particle equations
dh/dt = (hinf-h)/tauh : 1
hinf = 1/(1+exp((v-Vhinf)/khinf)) : 1
tauh = 1/(alphah + betah) : second
#betah = (-Abh*(v-Vbh))/(1-exp((v-Vbh)/kbh))/ms : Hz
alphah = Aah * kah/ exprel(-(v-Vah)/kah) /ms : Hz
betah = Abh * kbh/exprel((v - Vbh)/kbh)/ms : Hz
'''
paramKDR = dict()
# m particle
paramKDR["Aan"] = 0.02/mvolt # A paramenter of n alpha rate function
paramKDR["Abn"] = 0.002/mvolt # A parameter of n beta rate function
paramKDR["Van"] = 20*mvolt # V1/2 of alpha n
paramKDR["Vbn"] = 20*mvolt # V1/2 of beta n
paramKDR["kan"] = 9*mvolt # slope of alpha n
paramKDR["kbn"] = 9*mvolt # slope of beta n
eqIKDR = '''
#v : volt
# Current equation
IKDR = gKDR * n*(v-EK) : amp
# activation particle equations
dn/dt = alphan*(1-n)-betan*n : 1
alphan = Aan*kan/exprel(-(v - Van)/kan)/ms : Hz
betan = Abn*kbn/exprel((v-Vbn)/kbn)/ms : Hz
'''
# Hugenard slow K2 current
# Boltzman params
paramK2 = dict()
paramK2['Vlinf'] = -43*mvolt
paramK2['klinf'] = -17*mvolt
paramK2['Voinf'] = -58*mvolt
paramK2['koinf'] = 10.6*mvolt
# tau function params
paramK2['Vl1'] = 81 * mvolt
paramK2['kl1'] = 25.6 * mvolt
paramK2['Vl2'] = 132 * mvolt
paramK2['kl2'] = -18 * mvolt
paramK2['constl'] = 9.9 * msecond
paramK2['Vo1'] = 1329 * mvolt
paramK2['ko1'] = 200 * mvolt
paramK2['Vo2'] = 130 * mvolt
paramK2['ko2'] = -7.1 * mvolt
paramK2['consto'] = 120 * msecond
paramK2['constp'] = 8.9 * second
eqIK2 = '''
#v : volt # comment out if needed for CC experis
# Current equation
IKslow = 0.6*IK2a + 0.4*IK2b : amp
IK2a = l*o *gK2 *(v-EK): amp
IK2b = l*p*gK2 *(v-EK): amp
# activation particle equations
dl/dt = (linf-l)/(taul) : 1
linf = (1/(1+exp((v-Vlinf)/klinf)))**4 : 1
taul = 1/((exp((v - Vl1)/kl1) + exp((v + Vl2)/kl2)))*msecond + constl : second
# inactivation particle equations
do/dt = (oinf-o)/(tauo) : 1
oinf = 1/(1+exp((v-Voinf)/koinf)) : 1
tauo = 1/((exp((v - Vo1)/ko1) + exp((v + Vo2)/ko2)))*msecond + consto : second
# Slow inactivation equation, same exept for depolarised potentials
dp/dt = (oinf-p)/(taup+0.001*ms) : 1
taup = tauo*int(((v/mvolt) < -70)) + constp*int(((v/mvolt) > -70)) : second
'''
# General constants
A = 680 * umetre**2 # membrane area
Cconst = 1*ufarad*cm**-2 # cap per membrane area
Cm = A * Cconst # membrane Capacitance
ENa = 55*mvolt # Reversal potential Sodium
EK = -90*mvolt # K current reversal
# These were chosen just by eyeballing the resulting VC traces, do not take too seriously
densNaT = 4000*psiemens/pfarad # Density for NaT,t
# Density for KDR
densKDR = 1280*psiemens/pfarad
densK2 = 1000*psiemens/pfarad
# Density for H current roughly based on Liu et al from fig1
gNaT = densNaT*Cm # gmax for NaT etc
gKDR = densKDR*Cm # gmax for KDR
gK2 = densK2*Cm # gmax for fast inactivating K2
# set up the systen
eqsCC = '''
dv/dt = (Im + I)/Cm : volt
Im = (INaT + IKDR + IKslow) : amp
I : amp (shared) # stimulation current
'''
eqsCurrents = eqINaT + eqIKDR + eqIK2
eqsCCcombined = eqsCC + eqsCurrents
start_scope()
# collect all variables to one dict
allparam = {**paramNaT, **paramKDR, ** paramK2}
# make a neuron
SGNtest = NeuronGroup(1,
eqsCCcombined,
method='exponential_euler',
namespace=allparam)
# timestep/sampling rate
defaultclock.dt = 0.01*ms
# set voltage first, then set all gating variables to their ss value
SGNtest.v = -60*mvolt
# for INaT
SGNtest.m = '1/(1 + betam/alpham)'
SGNtest.h = 'hinf'
# # for K2
SGNtest.l = 'linf'
SGNtest.o = 'oinf'
SGNtest.p = 'oinf'
# for IKDR
SGNtest.n = '1/(1 + betan/alphan)'
statemon = StateMonitor(SGNtest, ['v', 'Im'], record=True)
groupnet = Network(SGNtest, statemon)
groupnet.run(1000*ms)
What you have aready tried
Plotted out all rate (tau) functions and steady state functions of the channels and inspected them (look fine as far as I can tell).
There should not be any definition gaps (most have a pattern of 1/(1 +exp(V-Vh/k)) or variants thereof and seem unsuspicious to me.
Set the particle position to steady state before starting the run to make the simulation start correctly.
Since this happens with all models using the parametrisation of the HH equations with tau whereas the channel models using the alpha and beta pattern are fine I think I must have gone wrong there.
Expected output (if relevant)
Nice neuron settling at RMP.
Actual output (if relevant)
DIV0 error or unexpectedly large value for all gating particle variables.
I and v of the simulation go toward infinity.