How to write a conductance based on the GHK equation

[luadam4c]

Hi guys,

Although I've been focusing on doing experiments the past few months, My P.I. is still very interested in using xolotl so I'm now back at it.

I'm trying to create custom conductances that match what I had in NEURON. However, I realized that cpplab expects all conductances to follow Ohm's law. If I want to model the T-type calcium current with the GHK current equation, how should I set g_bar and E?

Also, I'm a little confused by what is needed as method arguments and what is not. I've attached the MOD file and the draft HPP file I created by adapting from existing CaT files. Can you let me know how to modify things to make it work?

ITGHK_20190815.zip

Thanks, Adam

[sg-s]

@luadam4c I have implemented a CaT channel based on the 1998 Destexhe paper, which is available here (it's part of xolotl now)

I put it into a hodgkin-huxley model and it seems to work as expected, though you may want to tweak it a bit. In any case, I believe this answers your questions about how to integrate conductances where the current is nonlinear with the voltage, and also how to read out temperature, etc.

[sg-s]

You can also write the equation as is, and use the use_current = 1 flag in xolotl to integrate the voltage using currents -- this is written specifically for nonlinear conductances.

[luadam4c]

Thanks guys.

Alec, I saw that you created the Soplata files. On that note, do you have a working thalamocortical neuron model I can play with? I don't the what the appropriate g_bars should be and right now the default is set to zero.

Srinivas, thanks for your code. A few things I have yet to understand:

1. How are initial conditions for m and h set? Like what the NMODL INITIAL block does.
2. Where are variables like m, h, minf and hinf declared? Is there a list of variables that don't need to be declared in the custom conductance file?
3. How come I can't name my conductance with caps? When it's named "ITGHK" it produces a bunch of scope errors. But when it's named "iTGHK" those errors go away.
4. Not sure what you mean by "writing equation as is" when you set use_current = 1

Thanks, Adam

[luadam4c]

5. I tried using your Destexhe 1998 CaT.hpp When I add with no arguments

   x.soma.add('CaT')

   And then look at the properties

   x

   I see:

CaTDestexhe object (bf44496) with:

      pbar : NaN
      gbar : NaN
         h : NaN
         E : NaN
         m : NaN

How do I verify that default values have been set correctly?

[alec-hoyland]

I can answer a couple of these questions.

1. In general, V = -65, Ca = 0.05, and conductances default to steady-state (or should, I can't remember if we made that change yet already) https://xolotl.readthedocs.io/en/master/how-to/specify-initial-conditions/

2. All the properties and methods of conductances are described here https://xolotl.readthedocs.io/en/master/reference/c++/conductance/

Anything defined in the conductance header file doesn't need to be defined in the conductance specification. See here for a minimalist example.

3. This is a C++ quirk. I'm honestly not sure why -- but this is why we name our currents ACurrent and HCurrent instead of A and H. Something about namespaces.

Thanks guys.

Alec, I saw that you created the Soplata files. On that note, do you have a working thalamocortical neuron model I can play with? I don't the what the appropriate g_bars should be and right now the default is set to zero.

Srinivas, thanks for your code. A few things I have yet to understand:

1. How are initial conditions for m and h set? Like what the NMODL INITIAL block does.

2. Where are variables like m, h, minf and hinf declared? Is there a list of variables that don't need to be declared in the custom conductance file?

3. How come I can't name my conductance with caps? When it's named "ITGHK" it produces a bunch of scope errors. But when it's named "iTGHK" those errors go away.

4. Not sure what you mean by "writing equation as is" when you set use_current = 1

Thanks, Adam

[alec-hoyland]

@luadam4c Here are some thalamocortical models. https://github.com/alec-hoyland/xolotl-model-zoo

I would use model_howard or model_soplata, I think you'll have the best luck with those.

[sg-s]

@luadam4c

1. xolotl generally interprets NaN as "default value". So if you set anything to NaN, you are telling the C++ code to overwrite it with some default initial value. For m and h in a conductance, NaNs are overwritten with m_inf and h_inf. This is done when a conductance is connected to a compartment specifically, here

2. m, h, etc are decalred in conductance.hpp, specifically, here. You don't (and shouldn't) declare them again.

3. We've seen this -- we don't understand it, but it's some weird quirk of C++ namespaces.

4. Right now, the model you want to implement is such that the current

I = f(V, Ca)

where f is some nonlinear function of V. This is not in the "normal" Ohmic form

I = g*m*h(V-E)

and therefore you were understandably confused as to how xolotl would work with this. However, any nonlinear function of V can be rewritten in the Ohmic form:

I = g*V - g*E

where g may be some complex, time-varying function of V

By rewriting the ODE as such, we are effectively using the exponential Euler solvers that work for Ohmic conductances even here (This is not so wild -- even in the "Ohmic" linear conductances, m and h are complex nonlinear functions of V)

In the model I wrote, I rearranged the ODE.

What I mean is that you don't have to. You can write the ODE as it is in the Dexesthe paper and use the flag use_current = 1 and then xolotl will compute the currents and then use that to update the voltage.

[sg-s]

1. I tried using your Destexhe 1998 CaT.hpp When I add with no arguments

x.soma.add('CaT')

And then look at the properties

x

I see:

CaTDestexhe object (bf44496) with:

      pbar : NaN
      gbar : NaN
         h : NaN
         E : NaN
         m : NaN

How do I verify that default values have been set correctly?

The defaults probably aren't correct -- you should set them to sensible values. Note that in this model, gbar and E are not used (because the current equation does not depend on either term)

[luadam4c]

@luadam4c Here are some thalamocortical models. https://github.com/alec-hoyland/xolotl-model-zoo

I would use model_howard or model_soplata, I think you'll have the best luck with those.

Thanks, I'll definitely try them out!

[luadam4c]

• Variable declarations

  2. m, h, etc are decalred in conductance.hpp, specifically, here. You don't (and shouldn't) declare them again.

What about minf and hinf like in your code here?

    // integrate the gating variables
    minf = m_inf(V,Ca);
    m = minf + (m - minf)*exp(-dt/tau_m(V,Ca));

    hinf = h_inf(V,Ca);
    h = hinf + (h - hinf)*exp(-dt/tau_h(V,Ca));

When I changed that to mInf and hInf, it says it's not declared. So where are minf and hinf declared?

• The need for use_current

  4. Right now, the model you want to implement is such that the current

  I = f(V, Ca)

  where f is some nonlinear function of V. This is not in the "normal" Ohmic form

  I = g*m*h(V-E)

  and therefore you were understandably confused as to how xolotl would work with this. However, any nonlinear function of V can be rewritten in the Ohmic form:

  I = g*V - g*E

  where g may be some complex, time-varying function of V

  By rewriting the ODE as such, we are effectively using the exponential Euler solvers that work for Ohmic conductances even here (This is not so wild -- even in the "Ohmic" linear conductances, m and h are complex nonlinear functions of V)

  In the model I wrote, I rearranged the ODE.

  What I mean is that you don't have to. You can write the ODE as it is in the Dexesthe paper and use the flag use_current = 1 and then xolotl will compute the currents and then use that to update the voltage.

Wouldn't this part of the code cause it to update with currents already?

    container->i_Ca += g*V;

Or is the final integration really accessing just this

    g = pbar*m*m*h*G_by_V(V, Ca);

when use_current = 0?

• Also, I didn't understand why you had to write this method over again:

  double CaT::getCurrent(double V) {
  return V*G_by_V(V, container->Ca_prev);
  }

  You did not use it in integrate, so why define it? Also, even if you had to define it, would this be more correct?

  double CaT::getCurrent(double V) {
  return V*g;
  }

• Finally, what is the meaning of this line?

  gbar = gbar_next;

[luadam4c]

1. xolotl generally interprets NaN as "default value". So if you set anything to NaN, you are telling the C++ code to overwrite it with some default initial value. For m and h in a conductance, NaNs are overwritten with m_inf and h_inf. This is done when a conductance is connected to a compartment specifically, here

Then how do I verify whether the defaults are applied? For instance, if I have this in the hpp file:

if (isnan(pbar)) {pbar = .2e-3;}

and then add the conductance to a model

x.soma.add('CaT')

Shouldn't I get that default pbar value when I display the properties? But this is what I get:

pbar : NaN gbar : NaN h : NaN E : NaN m : NaN

[sg-s]

@luadam4c

verifying defaults:

they will be applied. Setting to NaN will guarantee this. To check, run it once, and NaN values will be overwritten.

[sg-s]

When I changed that to mInf and hInf, it says it's not declared. So where are minf and hinf declared?

minf is a variable that I use to store a value. m_inf is a function of voltage and calcium. You don't have to use minf, but if you do use it, you don't have to declare it (it's double). m_inf is a virtual function that you typically overwrite with your own function (but it must have the same function signature)

So the error you're getting is probably because you're changing the function signature of m_inf

[sg-s]

@luadam4c

Let's not worry about use_current for now because it's a needless complication at the moment.

To answer your question, yes, this will also work:

double CaT::getCurrent(double V) {
    return V*g;
}

When a compartment asks a conductance to integrate (under normal modes), it expects it to update the g variable (in the Ohmic case, this is g*m*h). The compartment then uses the updated g to update V -- it never computes the current! This follows the exponential Euler method in Dayan and Abbott (eq. 5.49)

The problem is with the conductance you are trying to model, we have an equation for the current that is in a non-Ohmic form. That's all right, because we can rewrite it so that it is in this form, and proceed as before.

[luadam4c]

When I changed that to mInf and hInf, it says it's not declared. So where are minf and hinf declared?

minf is a variable that I use to store a value. m_inf is a function of voltage and calcium. You don't have to use minf, but if you do use it, you don't have to declare it (it's double). m_inf is a virtual function that you typically overwrite with your own function (but it must have the same function signature)

So the error you're getting is probably because you're changing the function signature of m_inf

Nevermind, I found it. So minf is also declared in conductance.hpp at the end.

    // housekeeping, temp variables
    double minf = 0;
    double hinf = 1;
    double taum = 1;
    double tauh = 1;

Perhaps you can add these to the list here?

[luadam4c]

@luadam4c

verifying defaults:

they will be applied. Setting to NaN will guarantee this. To check, run it once, and NaN values will be overwritten.

I tried running it once. It still doesn't overwrite the NaNs

[luadam4c]

@luadam4c

Let's not worry about use_current for now because it's a needless complication at the moment.

To answer your question, yes, this will also work:

double CaT::getCurrent(double V) {
    return V*g;
}

When a compartment asks a conductance to integrate (under normal modes), it expects it to update the g variable (in the Ohmic case, this is g*m*h). The compartment then uses the updated g to update V -- it never computes the current! This follows the exponential Euler method in Dayan and Abbott (eq. 5.49)

The problem is with the conductance you are trying to model, we have an equation for the current that is in a non-Ohmic form. That's all right, because we can rewrite it so that it is in this form, and proceed as before.

So the getCurrent() method is just for the user to extract an individual component for plotting purposes?

So if V is updated by g, what gets updated by container->i_Ca? The calcium concentration? Which file is for the calcium decay mechanism then?

[alec-hoyland]

i_Ca is a property of the container of the conductance (i.e. the compartment) that keeps track of the instantaneous calcium current. i_Ca is used to compute the change in calcium concentration if you have a mechanism added to your model.

For example, to add the mechanism from Buchholtz et al. 1992, which is used in Destexhe et al. 1993 and others:

https://github.com/alec-hoyland/xolotl-model-zoo/blob/6edf91693f6dd30d2c117a9b245adb453294ae68/model_howard.m#L5

[sg-s]

@luadam4c

Yes, getCurrent is used only for plotting the currents in normal integration modes

Yes, V is updated only using g. container->i_Ca is used by a calcium mechanism (which is a different class, of type mechanism) to update the Ca in the compartment. If you don't have a calcium mechanism, your compartment won't change its Ca!

[luadam4c]

@luadam4c Here are some thalamocortical models. https://github.com/alec-hoyland/xolotl-model-zoo I would use model_howard or model_soplata, I think you'll have the best luck with those.

Thanks, I'll definitely try them out!

Alec,

model_soplata.m compiles as it is, but when I changed just the name of the compartment from comp to soma in model_soplata.m, as in the following:

x = xolotl('temperature', 36, 'temperature_ref', 22);

x.add('compartment', 'soma', 'Cm', 10, 'A', 1e-3, 'vol', 1e-6, 'Ca_out', 2e3);

x.soma.add('buchholtz/CalciumMech', 'Ca_in', 0.2, 'tau_Ca', 80, 'phi', 1);

x.soma.add('soplata/CaN', 'gbar', 0.25);
x.soma.add('destexhe/CaT', 'gbar', 1);
x.soma.add('soplata/Kd', 'gbar', 200);
x.soma.add('soplata/MCurrent', 'gbar', 7.5);
x.soma.add('soplata/NaV', 'gbar', 2000);
x.soma.add('Leak', 'gbar', 0.1);

I get the following error upon x.plot:

Error using X_a015ff8ae08dae31dad8eb38df052d34 gbars cannot be negative for any conductance

Error in xolotl/integrate (line 175) [results{1:n_outputs+1}] = f(arguments,self.I_ext',self.V_clamp');

Error in xolotl/plot (line 128) [V, Ca, ~, currents] = self.integrate;

Why is that?

[luadam4c]

Srinivas,

Starting with your code, I tried implementing our version of CaT.hpp with the temperature dependence included and added parameters for shifting the activation/inactivation curves. However, it is not converging when I voltage clamp. Can you see if I have any syntactic errors in this code?

// _  _ ____ _    ____ ___ _
//  \/  |  | |    |  |  |  |
// _/\_ |__| |___ |__|  |  |___
//
// T-type calcium conductance
//
// Adapted from:
//      Dendritic Low-Threshold Calcium Currents in Thalamic Relay Cells
//      Alain Destexhe, Mike Neubig, Daniel Ulrich and John Huguenard
//      https://www.jneurosci.org/content/18/10/3574.short

#ifndef CAT
#define CAT
#include "conductance.hpp"

// Define constants (compatible with NEURON)
//  https://www.neuron.yale.edu/phpBB/viewtopic.php?t=3783
const double FARADAY = 96485.309;
const double R = 8.3134;

// Inherit conductance class specifications
//  Note: Class names must not be all caps
class CaT_m3ha: public conductance {

private:

    // zzFF_by_RT = 4*F*F/(R*T)
    double zzFF_by_RT;

    // K_exp = exp(-ZFV/RT)
    double K_exp;

    // For temperature dependence
    double qm;
    double qh;
    double delta_temp;
    double phim;
    double phih;

public:

    // Note: the following are declared in conductance.hpp
    //      container, gbar, gbar_next, E, m, h, minf, hinf, taum, tauh, p, q, 
    //      dt, temperature, is_calcium

    // Maximum permeability in units of cm/s
    double pbar; 

    // Shifts in activation and inactivation curves
    double shiftm;
    double shifth;
    double slopem;
    double slopeh;

    // Specify parameters + Initial conditions
    CaT_m3ha(double g_, double E_, double m_, double h_, double pbar_, double shiftm_, double shifth_, double slopem_, double slopeh_)
    {
        // Read in arguments
        gbar = g_;
        E = E_;
        m = m_;
        h = h_;
        pbar = pbar_;
        shiftm = shiftm_;
        shifth = shifth_;
        slopem = slopem_;
        slopeh = slopeh_;

        // Note: gbar and E are unused in this model
        //  For safety, E must be set to zero if I = g(V-E) is accidentally used
        E = 0;

        // We can't set as NaN?
        gbar = 0;

        // Set defaults for arguments
        if (isnan(pbar)) { pbar = .2e-3; }  // maximum Ca++ permeability
        if (isnan(shiftm)) { shiftm = 1; }      // depolarizing shift of activation curve
        if (isnan(shifth)) { shifth = 1; }      // depolarizing shift of inactivation curve
        if (isnan(slopem)) { slopem = 1; }      // scaling factor for slope of activation curve
        if (isnan(slopeh)) { slopeh = 1; }      // scaling factor for slope of inactivation curve

        // Inaccessible parameters
        qm = 3.6;      // Q10 for activation, based on Coulter et al 1989
        qh = 2.8;      // Q10 for inactivation, based on Coulter et al 1989

        // specify exponents of m and h
        p = 2;
        q = 1;

        // Whether to allow this channel to be approximated
        approx_m = 0;
        approx_h = 0;

        // This is a calcium channel
        is_calcium = true;
    }

    // Method declarations
    string getClass(void);
    void connect(compartment *);
    void integrate(double, double);
    double GHK_by_Vp(double, double);
    double getCurrent(double);
    double m_inf(double, double);
    double h_inf(double, double);
    double tau_m(double, double);
    double tau_h(double, double);

};

string CaT_m3ha::getClass(){
    return "CaT_m3ha";
}

void CaT_m3ha::connect(compartment * pcomp_) {
    // Set the default gbar, m, h
    conductance::connect(pcomp_);

    // Compute the constant zzFF_by_RT = 4*F*F/(R*T)
    //  Note: units of [mol / umol] * [C / mol)] * [J/(V * mol)] / 
    //                      ([J/(K * mol)] * [K])
    //              = [C/(V * umol)]
    zzFF_by_RT = (1.0e6) * fast_pow(2, 2) * FARADAY * FARADAY / 
                    (R * (temperature + 273.16));

    // Compute temperature dependence
    delta_temp = ( temperature - 23 ) / 10;
    phim = pow(qm, delta_temp);
    phih = pow(qh, delta_temp);
}

double CaT_m3ha::getCurrent(double V) {
    return g * V;
}

void CaT_m3ha::integrate(double V, double Ca) {
    // Compute the new T channel maximal conductance
    // Note: Re-write the GHK equation in the form
    //          I = gV - gE
    //      where E = 0.
    //      units of [cm / s] * [C / (V * L)] * [10 mm / cm]
    //                  = [S * mm / L] * [L / 10^6 mm^3]
    //                  = [uS / mm^2]
    gbar = pbar * GHK_by_Vp(V, Ca) * (10.0);

    // Compute the activation gating variable at this time step
    minf = m_inf(V, Ca);
    taum = tau_m(V, Ca);
    m = minf + (m - minf) * exp(-dt / taum);

    // Compute the inactivation gating variable at this time step
    hinf = h_inf(V, Ca);
    tauh = tau_h(V, Ca);
    h = hinf + (h - hinf) * exp(-dt / tauh);

    // Compute the new T channel conductance
    // Note: this is used to integrate voltage
    g = gbar * m * m * h;

    // Add to the calcium current
    // Note: this updated the calcium concentration
    container->i_Ca += getCurrent(V);
}

double CaT_m3ha::GHK_by_Vp(double V, double Ca) {
    // Compute dimensionless constant K_exp = exp(-ZFV/RT)
    K_exp = exp(-V / (container->RT_by_nF));

    // Note: units of [C/(V * umol)] * [umol / L] = [C/(V * L)]
    return zzFF_by_RT * (Ca - container->Ca_out * K_exp) / (1 - K_exp);
}

double CaT_m3ha::m_inf(double V, double Ca)
{
    return 1.0 / ( 1.0 + exp( (V - shiftm + 57.0) / -6.2 * slopem) );
}

double CaT_m3ha::h_inf(double V, double Ca) 
{
    return 1.0 / (1.0 + exp( ( V - shifth + 81.0) / 4.0 * slopeh) ); 
}

double CaT_m3ha::tau_m(double V, double Ca) 
{
    return ( 0.612 + 1.0 / 
                ( exp( (V + 132.0 - shiftm) / -16.7) 
                + exp( (V + 16.8 - shiftm) / 18.2) ) ) 
            / (phim * slopem);

}

double CaT_m3ha::tau_h(double V, double Ca) 
{
    if ( V < -80 + shifth ) {
        return 1.0 * exp((V + 467 - shifth) / 66.6)
               / (phih * slopeh);
    } else {
        return ( 28 + 1.0 * exp( (V + 22 - shifth) / -10.5) ) 
               / (phih * slopeh);
    }
}

#endif

Again, this is based on this version of IT.mod:

INDEPENDENT {t FROM 0 TO 1 WITH 1 (ms)}

NEURON {
    SUFFIX IT
    USEION ca READ cai, cao WRITE ica
    GLOBAL qm, qh
    GLOBAL phim, phih
    RANGE pcabar, shiftm, shifth, slopem, slopeh
    RANGE minf, hinf, taum, tauh, ica
}

UNITS {
    (molar) = (1/liter)
    (mV)    = (millivolt)
    (mA)    = (milliamp)
    (mM)    = (millimolar)

    FARADAY = (faraday) (coulombs)
    R       = (k-mole) (joule/degC)
}

PARAMETER {
    : GLOBAL variables whose values are fixed
    qm      = 3.6   (1)     : Q10 for activation, based on Coulter et al 1989
    qh      = 2.8   (1)     : Q10 for inactivation, based on Coulter et al 1989

    : RANGE variables whose values are specified in hoc
    pcabar  = .2e-3 (cm/s)  : default maximum Ca++ permeability
    shiftm  = 1     (mV)    : depolarizing shift of activation curve
    shifth  = 1     (mV)    : depolarizing shift of inactivation curve
    slopem  = 1     (1)     : scaling factor for slope of activation curve
    slopeh  = 1     (1)     : scaling factor for slope of inactivation curve
}

ASSIGNED {
    : Variables that are assigned outside the mod file
    v               (mV)
    celsius         (degC)
    cai                (mM)    : intracellular [Ca++] (mM)
    cao                (mM)    : extracellular [Ca++] (mM)

    : GLOBAL variables that are assigned in the INITIAL block
    phim            (1)     : temperature adjustment to taum
    phih            (1)     : temperature adjustment to tauh

    : RANGE variables that are assigned in the INITIAL & DERIVATIVE blocks
    minf            (1)     : steady state value of activation gating variable
    hinf            (1)     : steady state value of inactivation gating variable
    taum            (ms)    : time constant for activation
    tauh            (ms)    : time constant for inactivation

    : RANGE variables that are assigned in the BREAKPOINT block
    ica             (mA/cm2): calcium current (mA/cm2) generated
}

STATE {
    m                (1)        : activation gating variable
    h                (1)        : inactivation gating variable
}

INITIAL {
    : Equations for taum & tauh were taken from Huguenard & McCormick, 1992 
    :   and were at 23 degC
    : Transformation to celsius using Q10 (qm & qh)

    : Calculate constants from temperature
    phim = qm ^ ( ( celsius - 23 (degC) ) / 10 (degC) )
    phih = qh ^ ( ( celsius - 23 (degC) ) / 10 (degC) )

    : Initialize state variables to steady state values
    evaluate_fct(v)            : calculate steady state values of state variables
    m = minf
    h = hinf
}

BREAKPOINT {
    : Update gating variables m & h from voltage v
    SOLVE castate METHOD cnexp

    : Calculate calcium current ica from m, h, cai, cao, v
    ica = pcabar * m*m*h * ghk(v, cai, cao)
}

DERIVATIVE castate {
    : Update minf, hinf, taum, tauh from v
    evaluate_fct(v)

    : Update m & h from minf, hinf, taum, tauh
    m' = (minf - m) / taum
    h' = (hinf - h) / tauh
}

PROCEDURE evaluate_fct(v (mV)) {
    : Update minf, hinf from v
    minf = 1.0 / ( 1 + exp( (v + 57 (mV) - shiftm) / (-6.2 (mV) * slopem) ) )
    hinf = 1.0 / ( 1 + exp( (v + 81 (mV) - shifth) / (4.0 (mV) * slopeh) ) )

    : Update taum, tauh from v
    taum = ( 0.612 (ms) + 1.0 (ms) / 
                                ( exp( (v + 132 (mV) - shiftm) / -16.7 (mV)) 
                                + exp( (v + 16.8 (mV) - shiftm) / 18.2 (mV)) ) ) 
            / (phim * slopem)
    if( v < -80 + shifth) {
        tauh = 1.0 (ms) * exp((v + 467 (mV) - shifth) / 66.6 (mV))
               / (phih * slopeh)
    } else {
        tauh = ( 28 (ms) + 1.0 (ms) * exp( (v + 22 (mV) - shifth) / -10.5 (mV)) ) 
               / (phih * slopeh)
    }
}

FUNCTION ghk(v (mV), ci (mM), co (mM)) (.001 coul/cm3) {
    LOCAL z, eci, eco

    z = (1e-3) * 2 * FARADAY * v / (R * (celsius + 273.15 (degC) ) )
                            : this is ZFV/RT, which is dimensionless 
                            : after applying conversion factor 1e-3 V/mV
    eco = co * efun(z)      : this has units of [mM] = [umol/cm3]
    eci = ci * efun(-z)     : this has units of [mM] = [umol/cm3]
    ghk = (1e-3) * 2 * FARADAY * (eci - eco)
                            : this has units of [mC/cm3]
                            : after applying conversion factor 1e-3 mmol/umol
}

FUNCTION efun(z) {
    if (fabs(z) < 1e-4) {
        efun = 1 - z/2      : 1st order Taylor approximation of z/(exp(z) - 1)
    } else {
        efun = z / (exp(z) - 1)
    }
}

FUNCTION nongat(v (mV), cai (mM), cao (mM)) (mA/cm2) {    : non gated current
    nongat = pcabar * ghk(v, cai, cao)
}

[sg-s]

1. The define guards name should match the class name: so on line 12, it should be "CaT_m3ha"

otherwise it looks good. as to why it doesn't converge with voltage clamp -- need to check

[sg-s]

I am able to voltage clamp the cell with the CaT conductance you're working with:

x = xolotl.examples.BurstingNeuron('prefix','liu');
x.AB.CaT.destroy()
x.AB.add('destexhe98/CaT')
x.AB.HCurrent.gbar=  20;

% now clamp
t = (1:x.t_end/x.sim_dt)*x.sim_dt;
x.V_clamp = 5*sin(t/100)' - 50;
I = x.integrate;
plot(I)

[sg-s]

I hope everything works for you, closing for now