Modeling leaky integrate-and-fire neuron including Hodgkin–Huxley-like currents

NMODL and the Channel Builder.
Post Reply
harrisonZ
Posts: 8
Joined: Thu Jan 23, 2020 3:45 pm

Modeling leaky integrate-and-fire neuron including Hodgkin–Huxley-like currents

Post by harrisonZ »

I have constructed a MOD to model leaky integrate-and-fire neuron including HH-like currents from the paper "Modeling the Effects of Transcranial Magnetic Stimulation on Cortical Circuits". The link to the paper is:
https://www.physiology.org/doi/full/10. ... b%3Dpubmed

The mechanism is governed by two equations: dv/dt for the membrane potential (v) and dtheta/dt for the firing threshold (theta).

Here is my code:

Code: Select all

TITLE esser_mech.mod 

NEURON {
    SUFFIX esser
    POINT_PROCESS esser_mech
    USEION na READ ena WRITE ina
    USEION k READ ek WRITE ik
    RANGE gNa_leak, gK_leak, theta_eq, tau_theta, tau_spike, tau_m, gspike, C, tspike, vrefrac, ifake
    NONSPECIFIC_CURRENT ispike
}
UNITS {
    (mV) = (millivolt)
    (mA) = (milliamp)
    (S) = (siemens)
}

PARAMETER {
    gNa_leak = 0.14 (S/cm2)
    gK_leak = 1.0 (S/cm2)
    C = 0.85 
    theta_eq = -53 (mV)  : resting threshold 
    tau_theta = 2 (ms) : threshold time constant
    tau_spike = 1.75 (ms) : time const during a spike
    tau_m = 15 (ms) : membrane time const
    tspike = 2 (ms) : duration of spike 
    vap = 30 (mV) : v of spike 
}
ASSIGNED {
    v (mV)
    ena (mV)
    ek (mV)

    ina (mA/cm2)
    ik (mA/cm2)
    ispike (mA/cm2)
    gspike (S/cm2)
    ifake (mA/cm2) : causes solver to execute the breakpoint block at the
:	same time as the calculation of the other currents
}

STATE {
    theta (mV)
}

BREAKPOINT {
    SOLVE states METHOD cnexp
    ina = gNa_leak*(v-ena)/tau_m     : Sodium-like leak current
    ik = gK_leak*(v-ek)/tau_m    : Potassium-like leak current
    ispike = gspike*(v-ek)/tau_spike 
    ifake = 0
}

INITIAL {
    net_send(0,1)
    theta = theta_eq
    gspike = 0
    ena = 30 (mV)
    ek = -90 (mV)
    v = -75.263 (mV)
}
DERIVATIVE states {
    theta' = (-(theta - theta_eq) + C*(v - theta_eq))/tau_theta  : threshold 
}

NET_RECEIVE (w) {
    if (flag == 1){
        WATCH (v > theta) 2
    }else if (flag == 2){
       net_event(t)
        gspike = 1 
        v = vap
        theta = vap 
        net_send(tspike,3)
    }else if (flag == 3){
        gspike = 0
        net_send(0,1)
    }
}
In the NET_RECEIVE block, I included the threshold detection using WATCH and spiking using net_event(t). gspike is 1 during an action potential and 0 when a cell is not spiking, and tau_spike is the time constant during a spike. I used net_send(tspike,3) to allow the AP to last tspike before stepping in to (flag==3).

My concerns are:
1. Is my implementation correct?
2. I know v is reserved for membrane potential, which is computed by numerical integration of the discretized cable equation. dv/dt equation in the paper (with a tau_m term) is not the same as the cable equation in NEURON. Would v be calculated correctly? Or do I need to define another variable for the membrane potential?

Thank you!
ted
Site Admin
Posts: 6287
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Modeling leaky integrate-and-fire neuron including Hodgkin–Huxley-like currents

Post by ted »

1. Is my implementation correct?
How will you discover whether it works properly? Did Esser et al. publish a figure that shows how one of their model cells responds to injected current or synaptic activation? Or is such a figure in their previous paper that they cited as the basis for this model? If not, you're cooked.

If you do find such a figure, your task has just begun. Esser et al. used point neuron models--model cells that have no spatial extent--so all of their currents and conductances are in "absolute" units, i.e. not current density or conductance density. Furthermore, density mechanisms can't have NET_RECEIVE blocks (unless that feature was added recently). So there's two reasons why you'll need to implement the mechanism as a POINT_PROCESS. And you'll have to decide what surface area your model cells should have, or neither synaptic mechanisms nor the "spike" mechanism will have appropritate effects on membrane potential.

I suspect that the WATCH statement in the NET_RECEIVE block

WATCH (v > theta) 2

is not going to do what you want. Will the threshold really vary, or should one regard a WATCH statement as being "set and forget" i.e. the threshold for spike detection becomes the value of theta at the instant that the WATCH statement is executed and does not change until the WATCH statement is executed again. I strongly suspect that the latter is the case, so you can't WATCH (v > theta) to detect occurrence of a spike. But a workaround is possible (see below).
2. I know v is reserved for membrane potential, which is computed by numerical integration of the discretized cable equation. dv/dt equation in the paper (with a tau_m term) is not the same as the cable equation in NEURON. Would v be calculated correctly? Or do I need to define another variable for the membrane potential?
The model cells are single compartment, so v is governed by a single ODE. Create a section, set nseg to 1, and you get a single ODE. Add a mechanism that generates some currents, and you get a single ODE in which the currents affect the time course of v.

The chief implementational problem is how to deal with the "varying threshold for spike detection." The answer is: declare a new ASSIGNED variable called whatever you like--maybe vaux (as in "auxiliary"). Calculate the value of vaux at each fadvance by inserting this

vaux = v - theta

as the last statement in the BREAKPOINT block. Then your WATCH statement can be

WATCH (vaux > 0) 2

and the "if (flag == 2) clause" needs to include the statement

vaux = vap


Final note: the HH mechanism involves 3 ODEs for the gating states (m, n, and h) for a total of 4 per cell (remember that v has its own ODE). The mechanism used by Esser et al. uses only one itself (for theta), for a total of 2 ODEs per cell. So in theory a network of Esser model cells might occupy half the space and run twice as fast as a network that uses HH-style mechanisms. But that ignores the contribution from synaptic mechanisms. A conductance based AMPAergic synaptic mechanism is described by a single ODE; conductance-based NMDA and GABA-Aergic mechanisms involve 2 ODEs. So even if one can take advantage of superposition (i.e. using a single synaptic mechanism instance to deal with multiple aferrent streams), and assuming only AMPAergic and GABAergic synapses, each model cell plus its synapses requires a total of 7 ODEs (for HH-style spike mechanism) vs. 5 (for Esser et al. spike mechanism). Bottom line: the Esser et al. mechanism isn't really going to save much space of computation time. I wonder who has actually reimplemented their model and used it in new work?
Post Reply