Thanks for your reply. Here is some code to recreate the error. The following code is based off the Mainen and Sejnowski code. In the code below, at the top there are two groups of conductances. One group contains the conductances from the original model with some h conductances added. Everything works with this set of conductances. However, when I alter the conductances to the other group. I get a segmentation error as shown in my previous posting above. However, if I comment out the two lines of code that implement the ih:
'soma { insert iH ghbar_iH = gh_soma } '
and
'insert iH ghbar_iH = gh_dend //inserted this h current' (within the dendritic_only sections}
The program works fine with the alternate set of currents.
Any thoughts about why adding in ih with an alternate set of conductances gives a segmentation error would be greatly appreciated.
Thanks again!
Corinne
Code: Select all
//This code is based on demofig1.hoc from Mainen and Sejnowski. This version does not have a gui interface, but instead automatically initializes and runs the code and prints the output voltage plot.
// Below are the values in the original Mainen code
// and some h current I have added to the model
gna_soma =20
gkv_soma =200
gca_soma =0.3
gkm_soma =0.1
gkca_soma =3.0
gna_dend =20
gkm_dend =0.1
gkca_dend =3
gca_dend =0.3
gkv_dend =0
gna_node =30000
gkv_axon =2000
gh_soma =.0001
gh_dend =.0001
//An altrnate set of conductances
/*
gna_soma =0
gkv_soma =0
gca_soma =1.73
gkm_soma =1
gkca_soma =1.63
gna_dend =0
gkm_dend =0
gkca_dend =0
gca_dend =0.307499
gkv_dend =0
gna_node =5000
gkv_axon =0
gh_soma =0
gh_dend =0
*/
//create objects for writing data
objref savv, savt
savv = new File()
savt = new File()
savv.wopen("TESTMainenV.dat")
savt.wopen("TESTMainenT.dat")
//-------------------
objref sh, st, axonal, dendritic, dendritic_only
// load_proc("nrnmainmenu")
create soma
access soma
tstop = 1000
steps_per_ms = 40
dt = 0.05
// --------------------------------------------------------------
// passive & active membrane
// --------------------------------------------------------------
ra = 150
global_ra = ra
rm = 30000
c_m = 0.75
cm_myelin = 0.04
g_pas_node = 0.02
v_init = -70
celsius = 37
Ek = -90
Ena = 60
// --------------------------------------------------------------
// Axon geometry
//
// Similar to Mainen et al (Neuron, 1995)
// --------------------------------------------------------------
n_axon_seg = 5
create soma,iseg,hill,myelin[2],node[2]
proc create_axon() {
create iseg,hill,myelin[n_axon_seg],node[n_axon_seg]
soma {
equiv_diam = sqrt(area(.5)/(4*PI))
// area = equiv_diam^2*4*PI
}
if (numarg()) equiv_diam = $1
iseg { // initial segment between hillock + myelin
L = 15
nseg = 5
diam = equiv_diam/10 // see Sloper and Powell 1982, Fig.71
}
hill {
L = 10
nseg = 5
diam(0:1) = 4*iseg.diam:iseg.diam
}
// construct myelinated axon with nodes of ranvier
for i=0,n_axon_seg-1 {
myelin[i] { // myelin element
nseg = 5
L = 100
diam = iseg.diam
}
node[i] { // nodes of Ranvier
nseg = 1
L = 1.0
diam = iseg.diam*.75 // nodes are thinner than axon
}
}
soma connect hill(0), 0.5
hill connect iseg(0), 1
iseg connect myelin[0](0), 1
myelin[0] connect node[0](0), 1
for i=0,n_axon_seg-2 {
node[i] connect myelin[i+1](0), 1
myelin[i+1] connect node[i+1](0), 1
}
}
// --------------------------------------------------------------
// Spines
// --------------------------------------------------------------
// Based on the "Folding factor" described in
// Jack et al (1989), Major et al (1994)
// note, this assumes active channels are present in spines
// at same density as dendrites
spine_dens = 1
// just using a simple spine density model due to lack of data on some
// neuron types.
spine_area = 0.83 // um^2 -- K Harris
proc add_spines() { local a
forsec $o1 {
a =0
for(x) a=a+area(x)
F = (L*spine_area*spine_dens + a)/a
L = L * F^(2/3)
for(x) diam(x) = diam(x) * F^(1/3)
}
}
proc init_cell() {
// passive
forall {
insert pas
Ra = ra
cm = c_m
g_pas = 1/rm
e_pas = v_init
}
// exceptions along the axon
forsec "myelin" cm = cm_myelin
forsec "node" g_pas = g_pas_node
// na+ channels
forall insert na
forsec dendritic gbar_na = gna_dend
forsec "myelin" gbar_na = gna_dend
hill.gbar_na = gna_node
iseg.gbar_na = gna_node
forsec "node" gbar_na = gna_node
// kv delayed rectifier channels
iseg { insert kv gbar_kv = gkv_axon }
hill { insert kv gbar_kv = gkv_axon }
soma { insert kv gbar_kv = gkv_soma }
//Ih current
soma { insert iH ghbar_iH = gh_soma }
// dendritic channels (SOMA IS ACTUALLY ONE OF THESE SECTIONS)
forsec dendritic {
insert km gbar_km = gkm_dend
insert kca gbar_kca = gkca_dend
insert ca gbar_ca = gca_dend
insert cad
}
//I added this for the dendrites only
forsec dendritic_only {
insert kv gbar_kv = gkv_dend //insert kv currents
insert iH ghbar_iH = gh_dend //inserted this h current
}
soma {
gbar_na = gna_soma
gbar_km = gkm_soma
gbar_kca = gkca_soma
gbar_ca = gca_soma
}
forall if(ismembrane("k_ion")) ek = Ek
forall if(ismembrane("na_ion")) {
ena = Ena
// seems to be necessary for 3d cells to shift Na kinetics -5 mV
vshift_na = -5
}
forall if(ismembrane("ca_ion")) {
eca = 140
ion_style("ca_ion",0,1,0,0,0)
vshift_ca = 0
}
}
//The following is doing most of the heavy lifting of the code.
// This is where they do the loading of the morphology.
// and inject the current.
proc load_3dcell() {
// $s1 filename
aspiny = 0
forall delete_section() // removes all sections from the accessed section list
xopen($s1)
access soma
dendritic = new SectionList()
// make sure no compartments exceed 50 uM length
forall {
diam_save = diam
n = L/50
nseg = n + 1
if (n3d() == 0) diam = diam_save
dendritic.append()
}
// show cell
// sh = new PlotShape()
// sh.size(-300,300,-300,300)
dendritic_only = new SectionList()
forsec dendritic dendritic_only.append()
soma dendritic_only.remove()
create_axon() //creates the axon but I dont think it is plotting
init_cell() //puts in the values for the conductances and such
if (!aspiny) add_spines(dendritic_only,spine_dens)
st=new IClamp(.5)
st.dur = 900
st.del = 5
}
//this has all the files with the morphology
proc fig1a() {
load_3dcell("cells/lcAS3.hoc")
st.amp = 0.05
}
proc fig1b() {
load_3dcell("cells/j7.hoc")
st.amp = 0.07
}
proc fig1c() {
load_3dcell("cells/j8.hoc")
st.amp = 0.1
}
proc fig1d() {
load_3dcell("cells/j4a.hoc")
st.amp = 0.2
}
//-------------------------------------------------------
//Initialize and Integrate
//graph
objref g
g= new Graph ()
g.size(0,tstop, -80,60)
g.addvar("soma.v(0.5)", 1,1,0.6,0.9, 2)
proc initialize() {
finitialize(v_init)
fcurrent()
}
proc integrate(){
while (t<tstop) {
fadvance()
g.plot(t)
savv.printf("%f\n",soma.v(.5))
savt.printf("%f\n",t)
}
g.flush()
}
/* ----this is the original file gui method----
xpanel("Figure 1")
xbutton("a. L3 Aspiny","fig1a()")
xbutton("b. L4 Stellate","fig1b()")
xbutton("c. L3 Pyramid","fig1c()")
xbutton("d. L5 Pyramid","fig1d()")
xpanel()
nrnmainmenu()
nrncontrolmenu()
newPlotV() //I think this is a prespecified function but isnt documented
---------------------------- */
//---INSTEAD i AM JUST GOING TO CALL A PROCESS
fig1a() //calls load_3dcell which loads aspiny
initialize()
integrate()
savv.close()
savt.close()