Two months ago, I implemented the two-compartment model (Ladenbauer 2019) in Neuron. Here is the link to the last post regarding the implementation: viewtopic.php?f=18&t=4370. The model uses its own variables vs and vd for the transmembrane potentials at soma and dendrite. Now, I am trying to build a small network that has multiple synaptic connections.
First, I tried Exp2Syn. I wanted to use the current from an Exp2Syn to feed into the cell soma and affect its voltage, but the voltage is calculated in the mod file, and the mod file doesn't accept currents, it accepts events with a weight. I used 'setpointer' to link the synaptic current to the Pointer in the mod file. It worked fine with a single synaptic connection. However, it did not work well for multiple synaptic connections since I have to know the number of synaptic connections and have enough Pointers in the mod file.
Then, I looked into the Izhikevich2007a mod (https://senselab.med.yale.edu/modeldb/s ... mod#tabs-2). It uses its own variable V for membrane potential calculation and includes its own synaptic mechanisms. I combined my neuronal dynamics and their synaptic mechanisms into one mod. However, I got a Segmentation Violation error when I was trying to make the single neuron fire (without any connections).
The error message was shown as below:
Code: Select all
oc>C:\nrn\bin\nrniv.exe: Segmentation violation
near line 1
{run()}
^
fadvance()
advance()
step()
continuerun(300)
and others
Code: Select all
COMMENT
Created by Zhihe Zhao in 2021
ENDCOMMENT
NEURON {
POINT_PROCESS twocomps
RANGE cs, vs, gi, gs, ge, vd, gd, cd, deltaT, delta, E, E1, E0, vt, vth, vr, iion, is, id, is0, id0, stds, stdd, freq, phase
RANGE isyn, gAMPA, gNMDA, gGABAA, gGABAB, tauAMPA, tauNMDA, tauGABAA, tauGABAB, eAMPA, eNMDA, eGABAA, eGABAB
RANGE deltat, t0, factor, vfactor1, vfactor2
}
UNITS {
(mV) = (millivolt)
(mA) = (milliamp)
(nA) = (nanoamp)
(S) = (siemens)
(F) = (farad)
(uS) = (microsiemens)
}
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)
E1 = 1 (mV/mm)
E0 = 0 (mV/mm)
vt = 10 (mV)
vth = 20 (mV)
vr = 0 (mV)
stds = 0 (mA)
stdd = 0 (mA)
is0 = 10E-9 (mA)
id0 = 5E-9 (mA)
freq=1 (Hz)
phase=0 (rad)
PI=3.14159265358979323846
tauAMPA = 5 (ms) : Receptor time constant, AMPA
tauNMDA = 150 (ms) : Receptor time constant, NMDA
tauGABAA = 6 (ms) : Receptor time constant, GABAA
tauGABAB = 150 (ms) : Receptor time constant, GABAB
eAMPA = 0 (mV)
eNMDA = 0 (mV)
eGABAA = -70 (mV)
eGABAB = -90 (mV)
vfactor1 = 80 (mV)
vfactor2 = 60 (mV)
}
ASSIGNED{
iion (mA)
refractory
is (mA)
id (mA)
i1 (mA)
i2 (mA)
E (mV/mm)
isyn (nA)
gAMPA (uS): AMPA conductance
gNMDA (uS): NMDA conductance
gGABAA (uS): GABAA conductance
gGABAB (uS): GABAB conductance
factor : Voltage factor used for calculating the current
t0 (ms): Previous time
deltat : Time step
}
: Initial conditions
INITIAL {
net_send(0,1)
vs = 0
vd = 0
refractory = 0
isyn = 0
gAMPA = 0
gNMDA = 0
gGABAA = 0
gGABAB = 0
t0 = t
deltat = 0
}
STATE {
vs (mV)
vd (mV)
}
: seed for white noise
PROCEDURE seed(x) {
set_seed(x)
}
BEFORE BREAKPOINT{
i1 = is0 + stds*normrand(0,1)*1E-9
i2 = id0 + stdd*normrand(0,1)*1E-9
}
BREAKPOINT {
deltat = t-t0 : Find time difference
SOLVE states METHOD cnexp
iion = gs*vs - ge*deltaT*exp((vs-vt)/deltaT)
is = i1
id = i2
E = E1*sin(2*PI*freq*t/1000+phase)+E0
: Receptor dynamics -- the correct form is gAMPA = gAMPA*exp(-delta/tauAMPA), but this is 30% slower and, in the end, not really any more physiologically realistic
gAMPA = gAMPA - deltat*gAMPA/tauAMPA : "Exponential" decays -- fast excitatory (AMPA)
gNMDA = gNMDA - deltat*gNMDA/tauNMDA : Slow excitatory (NMDA)
gGABAA = gGABAA - deltat*gGABAA/tauGABAA : Fast inhibitory (GABA_A)
gGABAB = gGABAB - deltat*gGABAB/tauGABAB : Slow inhibitory (GABA_B)
: Calculate current
factor = ((vs+vfactor1)/vfactor2)*((vs+vfactor1)/vfactor2)
isyn = gAMPA*(vs-eAMPA) + gNMDA*factor/(1+factor)*(vs-eNMDA) + gGABAA*(vs-eGABAA) + gGABAB*(vs-eGABAB)
t0=t : Reset last time so delta can be calculated in the next time step
}
: Calculate neuronal dynamics;
DERIVATIVE states {
vs' = 0.001*(gi*(vd-vs-delta*E)+is-iion+isyn*1E-6)/cs : eqn for Vs (soma)
vd' = 0.001*(gi*(vs-vd+delta*E)+id-gd*vd)/cd : eqn for Vd (dendrite)
}
: Input received
NET_RECEIVE (wAMPA(uS), wNMDA(uS), wGABAA(uS), wGABAB(uS)) {
INITIAL { wAMPA=wAMPA wNMDA=wNMDA wGABAA=wGABAA wGABAB=wGABAB} : Insanely stupid but required, otherwise reset to 0,
: Check if there is a spike
if (flag == 1) { : inputs integrated only when excitable
WATCH (vs>vth) 2 : Check if threshold has been crossed, and if so, set flag=2
}else if (flag == 2) {
net_event(t) : Send spike event
vs = vth
net_send(0,3)
gAMPA = 0 : Reset conductances -- not mentioned in Izhikevich's paper but necessary to stop things from exploding!
gNMDA = 0
gGABAA = 0
gGABAB = 0
}else if (flag == 3) { : ready to integrate again
vs = vr
net_send(0,1)
gAMPA = gAMPA + wAMPA
gNMDA = gNMDA + wNMDA
gGABAA = gGABAA + wGABAA
gGABAB = gGABAB + wGABAB
}
}