How to model spike and reset condition Ladenbauer 2019

Discussions of particular models.

Moderator: tom_morse

Post Reply
harrisonZ
Posts: 7
Joined: Thu Jan 23, 2020 3:45 pm

How to model spike and reset condition Ladenbauer 2019

Post by harrisonZ »

I am interested in writing the two-compartment neuron model (Weak electric fields promote resonance in neuronal spiking activity: Analytical results from two-compartment cell and network models,2019) in NEURON. The link: https://journals.plos.org/ploscompbiol/ ... bi.1006974. In the paper, they used python.

Image

The first and the second equations are for the dynamics of the somatic and the dendritic membrane potential, Vs and Vd. The reset condition is similar to the integrate-and-fire. The authors state that when Vs increases beyond VT (in equation 3), it diverges to infinity in finite time due to the exponential term.

I implemented the equations in the Mod file and modified the NET_RECEIVE block in IntFire.mod. I tried to reset Vs to Vr when it passes the threshold value. However, when Vs passes VT, I get an error message saying exp() out of range.

I suspect that when Vs goes beyond VT, NEURON cannot handle infinity value. Is there a way to work around this error?
harrisonZ
Posts: 7
Joined: Thu Jan 23, 2020 3:45 pm

Re: How to model spike and reset condition Ladenbauer 2019

Post by harrisonZ »

Code: Select all

NEURON {
    POINT_PROCESS twocomps
    RANGE cs, vs, gi, gs, ge, vd, gd, cd, deltaT, delta, E, E0, vt, vth, vr, iion, is,id
}

UNITS {
    (mV) = (millivolt)
    (mA) = (milliamp)
    (S) = (siemens)
}

PARAMETER {
    cs = 9.8868E-12 (F)
    cd = 2.8879E-11 (F)
    gi = 1.2126E-9 (S)
    gs = 2.4824E-10 (S)
    ge = 3.2956E-10 (S) 
    gd = 8.8163E-10 (S)
    deltaT = 1.5 (mV)
    delta = 0.324 (mm)  
    E = 1 (mV/mm)
    E0 = 0 (V/m)
    vt = 10 (mV)
    vth = 20 (mV)
    vr = 0 (mV)
    is = 0 (mA)
    id = 0 (mA)

}
ASSIGNED{
    iion (mA)
    refractory
    vaux (mV)
}

STATE {
    vs (mV)
    vd (mV)
}

BREAKPOINT {
    SOLVE states METHOD cnexp
    iion = gs*vs - ge*deltaT*exp((vs-vt)/deltaT)
    vaux = vs-vt
}

INITIAL {
    net_send(0,1)
    vs = 0 (mv)
    vd = 0 (mv)
    refractory = 0
}

DERIVATIVE states {
    vs' = (gi*(vd-vs-delta*E)+is-iion)/cs  : eqn for Vs (soma)
    vd' = (gi*(vs-vd+delta*E)+id-gd*vd)/cd : eqn for Vd (dendrite)

}

NET_RECEIVE (w) {
	if (flag == 1) { : inputs integrated only when excitable
        WATCH (vaux>-0.01) 2
    }else if (flag == 2) {
			vs = vth 
			net_send(0,3)
            net_event(t)
    }else if (flag == 3) { : ready to integrate again
		vs = vr
        net_send(0,1)
	}
}
Here I attached my code
ted
Site Admin
Posts: 5897
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: How to model spike and reset condition Ladenbauer 2019

Post by ted »

Before addressing the stability issue, here are some problems with the source code for twocomps that need to be fixed. These were discovered by checking the file with modlunit.

UNITS block: F is unknown. If you meant farad, insert the declaration
(F) = (farad)

INITIAL block: (mv) needs to be deleted from the statements
vs = ...
and
vd = ...

DERIVATIVE block: Assuming F means farad, the scale factor
(0.001)*
needs to be inserted into the right hand sides of the equations
vs' = ...
and
vd' = ...
ted
Site Admin
Posts: 5897
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: How to model spike and reset condition Ladenbauer 2019

Post by ted »

Now some comments about stability.

The authors concocted the mathematical form of the instantanously activating regenerative current (second term on the right hand side of their Eq. 3) simply to produce threshold and spiking behavior. Instantaneously activating regenerative currents can pose special challenges to the stability of numerical integration. Attempts to deal with this may result in a NEURON implementation that is stable at the cost of generating simulation results that differ quantitatively from the authors' original implementation. If you're lucky, the quantitative differences will not be so large as to cause qualitative differences which might destroy the utility of the NEURON implementation.

One modification that might fix the stability problem is to introduce a simple first order delay; this seems most likely to avoid destructive qualitative differences. Another is to substitute a different mathematical form for the voltage-dependence of the regenerative current; to me, this seems more problematical.

Let me know whether the following changes, which introduce a first order delay between "membrane potential" and the regenerative current, improve stability sufficiently. If problems persist, then I will propose a specific alternative form for the voltage dependence of m.

PARAMETER block: declare
taum = 0.1 (millisecond)
ASSIGNED block: declare
minf (1)
STATE block: declare
m (1)
BREAKPOINT block: change
iion = gs*vs - ge*deltaT*exp((vs-vt)/deltaT)
to
iion = gs*vs - ge*deltaT*m
INITIAL block: insert
rates()
m = minf
iion = gs*vs - ge*deltaT*m
DERIVATIVE block: insert
rates()
m' = (minf - m)/taum
Also add the following PROCEDURE block:

Code: Select all

PROCEDURE rates() {
  minf = exp((vs-vt)/deltaT)
}
harrisonZ
Posts: 7
Joined: Thu Jan 23, 2020 3:45 pm

Re: How to model spike and reset condition Ladenbauer 2019

Post by harrisonZ »

Hi Ted,

Thank you so much for the help! I managed to solve the issue. Vs resets to Vr (0 mV) after passing VT (10mV). Now I have some problems with setting Vs to Vth (20 mV) before resetting. I suspect there is something wrong with my NET_RECEIVE block. I set Vs to Vth after passing the VT but it would not do it correctly. I attached an image below:
Image

Code: Select all

NEURON {
    POINT_PROCESS twocomps
    RANGE cs, vs, gi, gs, ge, vd, gd, cd, deltaT, delta, E, E0, vt, vth, vr, iion, is,id
}

UNITS {
    (mV) = (millivolt)
    (mA) = (milliamp)
    (S) = (siemens)
    (F) = (farad)
}

PARAMETER {
    cs = 9.8868E-12 (F)
    cd = 2.8879E-11 (F)
    gi = 1.2126E-9 (S)
    gs = 2.4824E-10 (S)
    ge = 3.2956E-10 (S) 
    gd = 8.8163E-10 (S)
    deltaT = 1.5 (mV)
    delta = 0.324 (mm)  
    E = 1 (mV/mm)
    E0 = 0 (V/m)
    vt = 10 (mV)
    vth = 20 (mV)
    vr = 0 (mV)
    is = 0 (mA)
    id = 0 (mA)
    taum = 0.1 (millisecond)
}
ASSIGNED{
    iion (mA)
    refractory
    vaux (mV)
    minf (1)
}

STATE {
    vs (mV)
    vd (mV)
    m(1)
}

BREAKPOINT {
    SOLVE states METHOD cnexp
    iion = gs*vs - ge*deltaT*m
    vaux = vs-vt
}

INITIAL {
    rates()
    m = minf
    iion = gs*vs - ge*deltaT*m
    net_send(0,1)
    vs = 0 
    vd = 0 
    refractory = 0
}

DERIVATIVE states {
    rates()
    m'=(minf-m)/taum
    vs' = 0.001*(gi*(vd-vs-delta*E)+is-iion)/cs  : eqn for Vs (soma)
    vd' = 0.001*(gi*(vs-vd+delta*E)+id-gd*vd)/cd : eqn for Vd (dendrite)
}

PROCEDURE rates() {
  minf = exp((vs-vt)/deltaT)
}

NET_RECEIVE (w) {
	if (flag == 1) { : inputs integrated only when excitable
        WATCH (vs>vt) 2
    }else if (flag == 2) {
			vs = vth 
			net_send(0,3)
            net_event(t)
    }else if (flag == 3) { : ready to integrate again
		vs = vr
        net_send(0,1)
	}
}
harrisonZ
Posts: 7
Joined: Thu Jan 23, 2020 3:45 pm

Re: How to model spike and reset condition Ladenbauer 2019

Post by harrisonZ »

Hi Ted,

I found the issue in my code and fixed it. "WATCH (vs>vt) 2" should be "WATCH (vs>vth) 2". Thanks.
Post Reply