/* -------------------------------------------------------------- Multi-compartment simulations of neocortical neurons DEMO Z. F. Mainen and T. J. Sejnowski (1996) Influence of dendritic structure on firing pattern in model neocortical neurons. Nature 382: 363-366. Demo of Figure 1. updated: 11/1/96 author: Zach Mainen zach@salk.edu or zach@cshl.org Corrections to methods of Figure 1 misprinted in paper: 1. Time step = 25 usec, not 250 usec 2. I_Na rate functions are all shifted 5 mV negative m: alpha = 0.182(v+30)/(1-exp(-(v+30)/9)) beta = -0.124(v+30)/(1-exp((v+30)/9)) h: alpha = 0.024(v+45)/(1-exp(-(v+45)/5)) beta = -0.0091(v+70)/(1-exp((v+70)/5)) beta_inf = 1/(1+exp(v+60)/6.2) 3. I_Ca activation not inactivation is given first 4. I_Km rates alpha = 0.001 (v+30)/(1-exp(-(v+30)/9)) beta = 0.001 (v+30)/(1-exp((v+30)/9) 5. g_kca = 3 (pS um^-2) Correction to Figure 4: "transfer impedance" should be "transfer conductance" Z^-1 = I/V (uS) Minor revisions to GUI code by Ted Carnevale 2/2011 to use load_file rather than load_proc, and to control locations of RunControl panel, model selection panel, voltage axis graph, and PlotShape. -------------------------------------------------------------- */ load_file("nrngui.hoc") objref sh, st, axonal, dendritic, dendritic_only create soma access soma tstop = 1000 steps_per_ms = 40 dt = 0.025 // -------------------------------------------------------------- // 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 gna_dend = 20 gna_node = 30000 gna_soma = gna_dend gkv_axon = 2000 gkv_soma = 200 gca = .3 gkm = .1 gkca = 3 gca_soma = gca gkm_soma = gkm gkca_soma = gkca // -------------------------------------------------------------- // Axon geometry // // Similar to Mainen et al (Neuron, 1995) // -------------------------------------------------------------- n_axon_seg = 5 create soma,iseg,hill,myelin[2],node[2] // after http://www.neuron.yale.edu/phpbb/viewtopic.php?f=8&t=858 func secarea() { local x, sum sum = 0 for (x,0) sum += area(x) return sum } proc create_axon() { create iseg,hill,myelin[n_axon_seg],node[n_axon_seg] soma { equiv_diam = sqrt(secarea()/(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 } // dendritic channels forsec dendritic { insert km gbar_km = gkm insert kca gbar_kca = gkca insert ca gbar_ca = gca insert cad } 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 } } proc load_3dcell() { // $s1 filename aspiny = 0 forall delete_section() 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) */ // but control position of PlotShape window -- Ted Carnevale sh = new PlotShape(0) sh.size(-300,300,-299.522,299.522) sh.variable("v") {sh.view(-300, -299.522, 600, 599.043, 265, 369, 200.64, 200.32)} dendritic_only = new SectionList() forsec dendritic dendritic_only.append() soma dendritic_only.remove() create_axon() init_cell() if (!aspiny) add_spines(dendritic_only,spine_dens) st=new IClamp(.5) st.dur = 900 st.del = 5 } // GUI // run control panel { xpanel("RunControl") xvalue("Init","v_init", 1,"stdinit()", 1, 1 ) xbutton("Init & Run","run()") xbutton("Stop","stoprun=1") xvalue("Continue til","runStopAt", 1,"{continuerun(runStopAt) stoprun=1}", 1, 1 ) xvalue("Continue for","runStopIn", 1,"{continuerun(t + runStopIn) stoprun=1}", 1, 1 ) xbutton("Single Step","steprun()") xvalue("t","t", 2 ) xvalue("Tstop","tstop", 1,"tstop_changed()", 0, 1 ) xvalue("dt","dt", 1,"setdt()", 0, 1 ) xvalue("Points plotted/ms","steps_per_ms", 1,"setdt()", 0, 1 ) xvalue("Scrn update invl","screen_update_invl", 1,"", 0, 1 ) xvalue("Real Time","realtime", 0,"", 0, 1 ) xpanel(0,105) } // menu of models 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 } xpanel("Figure 1") xbutton("a. L3 Aspiny","fig1a()") xbutton("b. L4 Stellate","fig1b()") xbutton("c. L3 Pyramid","fig1c()") xbutton("d. L5 Pyramid","fig1d()") xpanel(159,481) // graph of v(0.5) vs. t objref g { g = new Graph(0) g.size(0,1000,-80,40) {g.view(0, -80, 1000, 120, 265, 105, 300.48, 200.32)} graphList[0].append(g) g.addexpr("v(.5)", 1, 1, 0.8, 0.9, 2) }