I am interested in writing the twocompartment neuron model (Weak electric fields promote resonance in neuronal spiking activity: Analytical results from twocompartment cell and network models,2019) in NEURON. The link: https://journals.plos.org/ploscompbiol/ ... bi.1006974. In the paper, they used python.
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 integrateandfire. 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?
How to model spike and reset condition Ladenbauer 2019
Moderator: tom_morse
Re: How to model spike and reset condition Ladenbauer 2019
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.8868E12 (F)
cd = 2.8879E11 (F)
gi = 1.2126E9 (S)
gs = 2.4824E10 (S)
ge = 3.2956E10 (S)
gd = 8.8163E10 (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((vsvt)/deltaT)
vaux = vsvt
}
INITIAL {
net_send(0,1)
vs = 0 (mv)
vd = 0 (mv)
refractory = 0
}
DERIVATIVE states {
vs' = (gi*(vdvsdelta*E)+isiion)/cs : eqn for Vs (soma)
vd' = (gi*(vsvd+delta*E)+idgd*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)
}
}

 Site Admin
 Posts: 5869
 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
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' = ...
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' = ...

 Site Admin
 Posts: 5869
 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
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 voltagedependence 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((vsvt)/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:
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 voltagedependence 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((vsvt)/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((vsvt)/deltaT)
}
Re: How to model spike and reset condition Ladenbauer 2019
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:
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:
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.8868E12 (F)
cd = 2.8879E11 (F)
gi = 1.2126E9 (S)
gs = 2.4824E10 (S)
ge = 3.2956E10 (S)
gd = 8.8163E10 (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 = vsvt
}
INITIAL {
rates()
m = minf
iion = gs*vs  ge*deltaT*m
net_send(0,1)
vs = 0
vd = 0
refractory = 0
}
DERIVATIVE states {
rates()
m'=(minfm)/taum
vs' = 0.001*(gi*(vdvsdelta*E)+isiion)/cs : eqn for Vs (soma)
vd' = 0.001*(gi*(vsvd+delta*E)+idgd*vd)/cd : eqn for Vd (dendrite)
}
PROCEDURE rates() {
minf = exp((vsvt)/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)
}
}
Re: How to model spike and reset condition Ladenbauer 2019
Hi Ted,
I found the issue in my code and fixed it. "WATCH (vs>vt) 2" should be "WATCH (vs>vth) 2". Thanks.
I found the issue in my code and fixed it. "WATCH (vs>vt) 2" should be "WATCH (vs>vth) 2". Thanks.