Adding Ih to the Mainen Sejnowski model

The basics of how to develop, test, and use models.
Post Reply
Corinne
Posts: 38
Joined: Wed Feb 09, 2011 7:13 pm

Adding Ih to the Mainen Sejnowski model

Post by Corinne »

Hi, I am attempting to add the Ih current available at http://senselab.med.yale.edu/modeldb/Sh ... 008\ih.mod to a variant of the Mainen and Sejnowski, patdemo code. I put ih.mod in the patdemo folder and recompiled everything by typing the nrnivmodl command.

I then inserted the following lines in the demofig1.hoc code:

Code: Select all

gh_soma=.001 
soma { insert iH  ghbar_iH = gh_soma }
This appears to work--when I run 'nrngui demofig1.hoc' the code runs. When I change the level of gh_soma by typing something like gh_soma=.0005 on the command line, the spiking behavior changes as expected.

However, I am trying to use a slight variation of the Mainen model. My version of the code works fine until I insert the same
soma { insert iH ghbar_iH = gh_soma } line of code, after which I get a segmentation error.

Code: Select all

NEURON -- Release 7.1 (359:7f113b76a94b) 2009-10-26
Duke, Yale, and the BlueBrain Project -- Copyright 1984-2008
See http://www.neuron.yale.edu/credits.html

loading membrane mechanisms from /usr2/ct/modelDB/patdemo/x86_64/.libs/libnrnmech.so
Additional mechanisms from files
 cad.mod ca.mod h.mod kca.mod km.mod kv.mod na.mod
	1 
	0 
	1 
	1 
/usr/local/nrn/x86_64/bin/nrniv: Segmentation violation
 in TESTMainen.hoc near line 272
 integrate()
            ^
        fadvance()
      integrate()

I would really appreciate any ideas on how to fix this.

Thanks!
Corinne
ted
Site Admin
Posts: 5810
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Adding Ih to the Mainen Sejnowski model

Post by ted »

Can't tell without seeing your new statement in the context of other relevant code.
Corinne
Posts: 38
Joined: Wed Feb 09, 2011 7:13 pm

Re: Adding Ih to the Mainen Sejnowski model

Post by Corinne »

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()

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

Re: Adding Ih to the Mainen Sejnowski model

Post by ted »

iH uses euler. NEVER use euler. Change its METHOD to cnexp and see if that takes care of the problem.
Corinne
Posts: 38
Joined: Wed Feb 09, 2011 7:13 pm

Re: Adding Ih to the Mainen Sejnowski model

Post by Corinne »

That did the trick. Thanks a bunch!
ted
Site Admin
Posts: 5810
Joined: Wed May 18, 2005 4:50 pm
Location: Yale University School of Medicine
Contact:

Re: Adding Ih to the Mainen Sejnowski model

Post by ted »

Thanks for using NEURON. And thanks for asking a question whose answer allowed me to stress something that many might otherwise be unaware of: SOLVE statements should never use METHOD euler. euler is prone to instability (the solution was blowing up abruptly!) and is not threadsafe. cnexp is the method to use for HH-style mechanisms (voltage-gated channels described by a system of ordinary differential equations).
Post Reply