HH Models Integration Errors - Support - Brian Simulator

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.

Full traceback of error (if relevant)