Fronto-parietal visuospatial WM model with HH cells (Edin et al 2007)

 Download zip file   Auto-launch 
Help downloading and running models
Accession:98017
1) J Cogn Neurosci: 3 structural mechanisms that had been hypothesized to underlie vsWM development during childhood were evaluated by simulating the model and comparing results to fMRI. It was concluded that inter-regional synaptic connection strength cause vsWM development. 2) J Integr Neurosci: Given the importance of fronto-parietal connections, we tested whether connection asymmetry affected resistance to distraction. We drew the conclusion that stronger frontal connections are beneficial. By comparing model results to EEG, we concluded that the brain indeed has stronger frontal-to-parietal connections than vice versa.
Reference:
1 . Edin F, Macoveanu J, Olesen P, Tegnér J, Klingberg T (2007) Stronger synaptic connectivity as a mechanism behind development of working memory-related brain activity during childhood. J Cogn Neurosci 19:750-60 [PubMed]
2 . Edin F, Klingberg T, Stödberg T, Tegnér J (2007) Fronto-parietal connection asymmetry regulates working memory distractibility. J Integr Neurosci 6:567-96 [PubMed]
Citations  Citation Browser
Model Information (Click on a link to find other models with that property)
Model Type: Realistic Network;
Brain Region(s)/Organism: Neocortex;
Cell Type(s): Neocortex U1 L2/6 pyramidal intratelencephalic GLU cell; Abstract Wang-Buzsaki neuron;
Channel(s):
Gap Junctions: Gap junctions;
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: NEURON;
Model Concept(s): Working memory; Attractor Neural Network;
Implementer(s):
Search NeuronDB for information about:  Neocortex U1 L2/6 pyramidal intratelencephalic GLU cell;
/*
* This template defines a laboratory with a cell, and the functions allows
* the user to calculate single cell properties such as fI-curves.
*
* Author: Fredrik Edin, 2003.
* Address: freedin@nada.kth.se
*
*/

load_file( "ECell.hoc" )
load_file( "ICell.hoc" )

begintemplate LabCell

    /* Objects */
    objref electrode  // The electrode to the cell (not used now)
    objref cell      // The cell in the laboratory, 1 or 3 compartments
    objref ns         // Netstim, a stimulating spike train
    objref nss[1]     // Vector of stimulating spike trains
    objref exinvec    // Vector with synapses for every spike train (st) in nss
    objref intervec   // Vector with intervals between spikes for st in nss
    objref nspikevec  // Vector with no of spikes for st in nss
    objref wvec       // Vector with synaptic weights st in nss
    objref shiftvec   // Vector with phase shifts between st in nss
    objref noisevec   // Vector with noise level for st in nss ( between 0-1 )
    objref samevec    // Vector where element e at index i says that indices i 
                      // and e should have the same source
    objref cvode      // The variable time step integrator used
    objref gvec       // Vector with GABA currents over time
    objref avec       // Vector with AMPA currents over time
    objref nvec       // Vector with NMDA currents over time
    objref Xgvec      // Vector with external GABA currents over time
    objref Xavec      // Vector with external AMPA currents over time
    objref Xnvec      // Vector with external NMDA currents over time
    objref tvec[6]    // Five vectors for gvec, avec, nvec, Xgvec, Xavec, Xnvec
    objref cvec       // The sum of gvec, avec, nvec, Xgvec, Xavec, Xnvec
    objref vPlot      // A plot (for time-dependent variables)
    objref vPlot2     // Another plot (for time-dependent variables)
    objref xyPlot     // A third plot (any graph)
    objref Ivec       // Vector with applied currents
    objref fIvec      // Vector with frequencies at currents in Ivec
    objref Ifvec      // Vector with external frequencies at frequencies fIvec
    objref xx         // Vector with abscissa values for xyPlot
    objref yy         // Vector with ordinate values for xyPlot
    objref spikevec[1]// Spike vector vector
    objref tmp
    objref tmp2
    objref sn
    strdef str        // Temporary string
    strdef pl         // Temporary string (plot or not?)
    strdef pl1        // Temporary string (possibility to plot or not?)
    objref file       // Temporary file object
    
    /* Public functions */
    public fICurve    // Calculate the fI-curve of the cell
    public effFICurve // Calculate the fI-curve of the cell
    public PSP        // Stimulate the cell with a PSP
    public PSPtrain   // Stimulate the cell with a train of PSPs
    public PSPtrain2  // Stimulate the cell with a train of PSPs
    public PSPtrains  // Stimulate the cell with several trains of PSPs
    public oneSpike   // Evoke a spike
    public iClamp     // Apply an I-clamp current
    public PSC        // Calculate the total charge in a PSC
    public plotXY     // Plot x and y values contained in vector arguments
    public resetWvec
    
    /* Public objects or variables */
    public eCell      // The cell in the laboratory
    public cvode      // The numerical time step integrator
    public charge     // Total charge of an EPSC
    public gabaChrg   // Charge through an GABA channel during a PSC
    public ampaChrg   // Charge through an AMPA channel during a PSC
    public nmdaChrg   // Charge through an NMDA channel during a PSC
    public XgabaChrg  // Charge through an external GABA channel during a PSC
    public XampaChrg  // Charge through an external AMPA channel during a PSC
    public XnmdaChrg  // Charge through an external NMDA channel during a PSC
    public gvec       // Vector with GABA currents over time
    public avec       // Vector with AMPA currents over time
    public nvec       // Vector with NMDA currents over time
    public Xgvec      // Vector with external GABA currents over time
    public Xavec      // Vector with external AMPA currents over time
    public Xnvec      // Vector with external NMDA currents over time
    public cvec       // The sum of gvec, avec, nvec, Xavec, Xnvec
    public tvec       // Five vectors for gvec, avec, nvec, Xavec, Xnvec
    public vPlot      // A plot (for time-dependent variables)
    public vPlot2     // Another plot (for time-dependent variables)
    public xyPlot     // Another plot (for any variables)
    public Ivec, fIvec, Ifvec // Vectors for fI-curve
    public nss        // A vector with netstim objects
    public spikevec   // Vector recording spikes
    public exinvec    // Vector needed to specify parameters for synaptic  
    public intervec   // input into the cell
    public nspikevec  
    public wvec    
    public noisevec 
    public samevec 
    public shiftvec
     
    /* Creates an ECellTest object
    *
    *(Arg 1, var_dt:  0 or no argument = use cvode
    *                               >0 = fix dt)
    *(Arg 2, pl      : plot = plot, Other arguments or no arguments = no plot)
    *(Arg 3, stats   : Statistics printed if there are more than two args.) 
    *(Arg 4, celltype: 1 or nothing = E-Cell, 2 = I-Cell)
    *(Arg 5, E       : Absolute tolerance of the CVode)
    */
    proc init() { local i, n1, n2
	
	/* Cell parameters */
	if( numarg() > 3 ) {
	    celltype = $4
	} else {
	    celltype = 1
	}
	if( celltype == 1) {
	    cell = new ECell()
	} else {
	    cell = new ICell()
	}
	
	if( numarg() > 0 ) {
	    var_dt = $1
	} else {
	    var_dt = 0 
	}
	cvode = new CVode(1)
	cvode.active(var_dt <= 0)
	if( var_dt > 0 ) {
	    dt = var_dt
	}
	print "Cvode active: ", cvode.active()
	if( numarg() > 1 ) {
	    pl1 = $s2
	} else {
	    pl1 = "no plot"
	}
	if( numarg() > 2 ) {
	    printStat = 1
	} else {
	    printStat = 0
	}
	
	if( strcmp( pl1, "plot" ) == 0 ) {
	    plotPoss = 1
	} else {
	    plotPoss = 0
	}
	
	if( numarg() > 4 && cvode.active() ) {
	    E = $5
	    cvode.atol( E )
	}
	
	if( plotPoss ) {
	    
	    /* Voltage graphs */
	    vPlot = new Graph(0)
	    vPlot.view( 100, -80, 200, 20, 0, 20, 600, 230 )
	    vPlot.erase()
	    vPlot.label(1,0) /* Set position for next label */
	    vPlot2 = new Graph(0)
	    vPlot2.view( 100, -80, 200, 20, 0, 320, 600, 230 )
	    vPlot2.erase()
	    vPlot2.label(1,0) /* Set position for next label */
	    
	    /* Voltages */
	    str = "cell.getV()"
	    vPlot.addexpr( str, 1, 1 ) /* Black */
	    
	    /* Intracellular calcium */
	    str =  "cell.getCai(0)"
	    //vPlot2.addexpr( str, 1, 1 ) /* Black */
	    str = "cell.getCai(2)"
	    //vPlot2.addexpr( str, 2, 1 ) /* Red */
	    str =  "cell.getGCan()"
	    //vPlot2.addexpr( str, 3, 2 ) /* Black */
	    str =  "cell.getmCan()"
	    //vPlot2.addexpr( str, 4, 1 ) /* Red */
	    
	    /* Synaptic states */
	    str = "cell.syn[2].s"
	    //vPlot2.addexpr( str, 2, 1 ) /* Red */
	    str = "cell.syn[4].s"
	    //vPlot2.addexpr( str, 4, 1 ) /* Green */
	    
	    /* Synaptic weights */
	    n1 = 3
	    n2 = 2
	    for i = n1, n2 {
		sprint( str, "cell.getWsyn(%d,0,1)+70", i )
		//vPlot.addexpr( str, 4, 1 ) /* Green */
	    }
	    
	    /* Channel variables */
	    /* Soma */
	    /* Na g */
	    sprint( str, "cell.getINa()" )
	    //vPlot2.addexpr( str, 3, 2 ) /* Blue */
	    /* K i */
	    sprint( str, "cell.getIK()" )
	    //vPlot2.addexpr( str, 2, 2 ) /* Red */
	    /* Synaptic currents */
	    sprint( str, "cell.getIsyn(0,0)" )
	    vPlot2.addexpr( str, 1, 2 )
	    sprint( str, "cell.getIsyn(3,2)" )
	    vPlot2.addexpr( str, 2, 2 )
	    sprint( str, "cell.getIsyn(2,1)" )
	    vPlot2.addexpr( str, 3, 2 )
	    sprint( str, "cell.getIsyn(2,2)" )
	    vPlot2.addexpr( str, 4, 2 )
	    	    
	    xyPlot = new Graph(0)
	    xyPlot.view( 0, 0, 100, 100, 640, 270, 300, 230)
	    xyPlot.erase()
	}
	print "ECellTest.init"
    }
    
    proc reGraph() {
	vPlot.begin()
	vPlot2.begin()
    }
    
    /* produces an fI-curve 
    *
    * Arg 1, Lo   : Lower stimulus intensity in uA/cm2
    * Arg 2, Hi   : Higher stimulus intensity in uA/cm2
    * Arg 3, n    : No of data points. 
    *(Arg 4, dur  : Stimulus length
    *(Arg 5, nCell: number of afferent cells)
    *(Arg 6, exin : 0 = Gaba, 1 = Ampa, 2 = Nmda)
    *(Arg 6, w    : Weight in mS/cm2)
    */
    proc fICurve() { local dur, maxi, i, lo, hi, n, imu, ext
	
	lo = $1
	hi = $2
	n = $3
	if( numarg() > 3 ) {
	    dur = $4
	} else {
	    dur = 1000
	}
	if( numarg() > 4 ) {
	    nCell = $5
	    exin = $6
	    w = $7
	} else {
	    nCell = 1
	    exin = 1
	    w = 0.016
	}
	
	fIvec = new Vector( n )
	Ifvec = new Vector( n )
	Ivec = new Vector( n )
	Ivec.indgen( lo, hi, (hi-lo)/(n-1) )

	/* Runs simulations of duration dur at different currents */
	for i = 0, n - 1 {
	    print "Iter #", i+1, "a"
	    imu = Ivec.x[i] /* Current in uA/cm2 */
	    iClamp( imu, dur, "plot" )
	    //PSPtrain2( imu, 1000, dur, "plot" )
	    fIvec.x[i] = nspk / ( dur / 1000 )
	    /*print "Iter #", i+1, "b"
	    if( fIvec.x[i] == 0 ) {
		Ifvec.x[i] = 0
	    } else {
		PSPtrain( exin, w, 1000/(fIvec.x[i]*nCell), fIvec.x[i]*nCell*(dur/1000), 1, "plot" )
		Ifvec.x[i] = -curr
	    }*/
	}
		
	if( plotPoss ) {
	    plotXY( Ivec, fIvec )
	    plotXY( Ifvec, fIvec )
	}
	file = new File()
	file.wopen("I-ECell1Comp.vec")
	Ivec.vwrite( file )
	file.close()
	file.wopen("f-ECell1Comp.vec")
	fIvec.vwrite( file )
	file.close()
	file.wopen("If-ECell1Comp.vec")
	//Ifvec.vwrite( file )
	file.close()
    }
    
    /* produces an effective fI-curve (current source is synaptic activity
    * and fI-curve is calculated with background synaptic activity) 
    *
    * Arg 1, Lo   : Lower stimulus intensity in Hz ( x 1024 for 1024 cells )
    * Arg 2, Hi   : Higher stimulus intensity in Hz ( x 1024 for 1024 cells )
    * Arg 3, n    : No of data points. 
    * Arg 4, bk   : Background strength (weight, f = 1000 always)
    * Arg 5, w    : Stimulus weight
    * Arg 6, N    : No of cells
    *(Arg 7, dur  : Stimulus length
    *(Arg 8, exin : 0 = Gaba, 1 = Ampa, 2 = Nmda(default))
    *(Arg 9, wI   : GABA weight
    */
    proc effFICurve() { local dur, maxi, i, lo, hi, n, imu, ext, del
	
	lo = $1
	hi = $2
	n = $3
	bk = $4
	w = $5
	N = $6
	if( numarg() > 6 ) {
	    dur = $7
	} else {
	    dur = 1000
	}
	if( numarg() > 7 ) {
	    exin = $8
	} else {
	    exin = 2
	}
	if( numarg() > 8 ) {
	    wI = $9
	} else {
	    wI = -1
	}
	
	if( celltype == 1 && exin > 0 ) {
	    m = 3 
	} else {
	    m = 2
	}
	m = m + (wI>0)
	
	exinvec = new Vector(m)
	intervec = new Vector(m)
	nspikevec = new Vector(m)  
	wvec = new Vector(m)
	shiftvec = new Vector(m-1) 
	noisevec = new Vector(m) 
	samevec = new Vector(m,-1)
	start = 150
	del = 1000
	
	for i = 0, m-2 { 
	    noisevec.x[i] = 1
	    shiftvec.x[i] = 0
	    wvec.x[i] = w
	}
	noisevec.x[m-1] = 1
	wvec.x[m-1] = w
	if( wI>0 ) {
	    wvec.x[m-1] = wI
	    intervec.x[m-1] = 1000/3*256
	    nspikevec.x[m-1] = (dur+del)/intervec.x[m-1]
	    exinvec.x[m-1] = 0
	    samevec.x[m-1] = -1
	}
	
	intervec.x[0] = 1000/1000
	nspikevec.x[0] = dur+del
	wvec.x[0] = bk
	
	if( celltype == 1 && exin > 0 ) {
	    exinvec.x[0] = 16
	    exinvec.x[1] = 15 + exin
	    exinvec.x[2] = 9 + exin
	    samevec.x[2] = 1
	} else if( celltype == 1 ) {
	    exinvec.x[0] = 16
	    exinvec.x[1] = 3 + exin
	} else {
	    exinvec.x[0] = 4
	    exinvec.x[1] = 3 + exin
	}
	
	fIvec = new Vector( n )
	Ivec = new Vector( n )
	Ivec.indgen( lo, hi, (hi-lo)/(n-1) )

	/* Runs simulations of duration dur at different currents */
	for i = 0, n - 1 {
	    print "---------------------"
	    print " "	    
	    print "Iter #", i+1, "a"
	    //imu = Ivec.x[i] /* Current in uA/cm2 */
	    //iClamp( imu, dur, "plot" )
	    //PSPtrain2( imu, 10000, dur, "plot" )
	    if( Ivec.x[i] == 0 ) {
		nspikevec.x[1] = 0
		intervec.x[1] = 1
	    } else {
		intervec.x[1] = 1000 / (N*Ivec.x[i])
		nspikevec.x[1] = (del+dur)/intervec.x[1]
	    }
	    fIvec.x[i] = PSPtrains( start, del, exinvec, intervec, nspikevec, shiftvec, noisevec, wvec, samevec )
	    print "Frequency: ", Ivec.x[i], " Weight: ", w
	    print "Firing rate: ", fIvec.x[i], " Hz"
	    //fIvec.x[i] = nspk / ( dur / 1000 )
	    if( plotPoss && i > 0 ) {
		plotXY( Ivec.c(0,i), fIvec.c(0,i) )
		//tmp = new Vector(2)
		//tmp.indgen(fIvec.x[0], Ivec.x[i])
		tmp2 = new Vector(2)
		tmp2.indgen(Ivec.x[0], Ivec.x[i])
		plotXY( tmp2, tmp2 )
	    }
	}
		
	file = new File()
	file.wopen("I-ECell1Comp.vec")
	Ivec.vwrite( file )
	file.close()
	file.wopen("f-ECell1Comp.vec")
	fIvec.vwrite( file )
	file.close()
    }
    
    /* Produce a PSP
    *
    * Arg 1, exin : 0 = GABA synapse
    *               1 = AMPA synapse
    *               2 = NMDA synapse
    * Arg 2, w    : The synaptic weight in mS/cm2
    *(Arg 3, ext  :  -1 = both local and external
    *                 0 = local
    *                 1 = external
    */
    proc PSP() { local exin, w, i
	exin = $1
	w = $2
	if( numarg() > 2 ) {
	    ext = $3
	} else {
	    ext = -1
	}
	cell.soma ns = new MyNetStim(.5)
	ns.interval = 0
	ns.number = 1
	ns.start = 150
	ns.noise = 0
	if( ext ) {
	    cell.pre_list[exin+3].append( new NetCon( ns, cell.syn[exin+3], 0, delay, 0 ) )
	    cell.setWsyn( exin+3, w )
	}
	if( ext < 1 ) {
	    cell.pre_list[exin].append( new NetCon( ns, cell.syn[exin], 0, delay, 0 ) )
	    cell.setWsyn( exin, w )
	}
	//cell.connect_pre( ns, 0, exin, 1 )
	cell.activate_syn( 1, exin )
	tStop = 400
	if( plotPoss ) {
	    vPlot.size( 0, tStop, -80, 20 )
	    vPlot2.size( 150, 155, 0, 0.001 )
	    spikevec = new Vector()
	    cell.pre_list[exin+(ext>0)*3].object(0).record( spikevec )
	    run("plot")
	    for i = 0, spikevec.size()-1 {
		vPlot2.mark( spikevec.x[i], 0.0006, "o", 2 )
	    }
	} else {
	    run("no plot")
	}
    }
    
    /* Produce a PSC. Charge flowing through different synapses is recorded
    *
    * Arg 1, gaba : gaba conductance in mS/cm2
    * Arg 2, nmda : nmda conductance in mS/cm2
    * Arg 3, ampa : ampa conductance in mS/cm2
    * Arg 4, Xgaba: external ampa conductance in mS/cm2
    * Arg 5, Xampa: external ampa conductance in mS/cm2
    * Arg 6, Xnmda: external nmda conductance in mS/cm2
    *(Arg 7, pl   : plot = plot, other arguments or no argument = don't plot. )
    *(Arg 8, disp : 1 = display data on screen, 0 = do not display)
    * Return charge : Total PSC charge in uC/cm2
    */
    xopen( "../../NeuronKlasser/integral.hoc" )
    func PSC() { local gaba, nmda, ampa, Xampa, Xnmda, i, low, disp
	
	gaba = $1
	ampa = $2
	nmda = $3
	Xgaba = $4
	Xampa = $5
	Xnmda = $6
	nCell = $7
	tStop = 1500 //FIX ME hitta lagom längd
	if( numarg() > 7 ) {
	    pl = $s8
	} else {
	    pl = "no plot"
	}
	if( numarg() > 7 ) {
	    disp = $8
	} else {
	    disp = 1
	}
	if( disp ) { 
	    print " "
	    print "ECellTest.PSC"
	}	
	
	cell.soma ns = new NetStim(.5)
	ns.interval = 0
	ns.number = 1
	ns.start = 150
	ns.noise = 0
	cell.connect_pre( ns, 0, 0, 1 )
	cell.connect_pre( ns, 0, 1, 1 )
	cell.connect_pre( ns, 0, 2, 1 )
	cell.pre_list[3].append( new NetCon( ns, cell.syn[3], 0, delay, 0 ) )
	cell.pre_list[4].append( new NetCon( ns, cell.syn[4], 0, delay, 0 ) )
	cell.pre_list[5].append( new NetCon( ns, cell.syn[5], 0, delay, 0 ) )
	cell.setWsyn( 0, gaba )
	cell.setWsyn( 1, ampa )
	cell.setWsyn( 2, nmda )
	cell.setWsyn( 0, Xgaba, file, 0, 1 )
	cell.setWsyn( 1, Xampa, file, 0, 1 )
	cell.setWsyn( 2, Xnmda, file, 0, 1 )
	cell.activate_syn( 1 )
	if( cvode.active() ) {
	    for i = 0, 5 {
		tvec[i] = new Vector(10000)
	    }
	    cvec = new Vector(10000)
	    gvec = new Vector(10000)
	    avec = new Vector(10000)
	    nvec = new Vector(10000)
	    Xgvec = new Vector(10000)
	    Xavec = new Vector(10000)
	    Xnvec = new Vector(10000)
	    cvode.record( &cell.syn[0].i, gvec, tvec[0] )
	    cvode.record( &cell.syn[1].i, avec, tvec[1] )
	    cvode.record( &cell.syn[2].i, nvec, tvec[2] )
	    cvode.record( &cell.syn[3].i, Xavec, tvec[3] )
	    cvode.record( &cell.syn[4].i, Xnvec, tvec[4] )
	    cvode.record( &cell.syn[5].i, Xnvec, tvec[5] )
	} else {
	    cvec = new Vector( tStop / dt )
	    gvec = new Vector( tStop / dt )
	    avec = new Vector( tStop / dt )
	    nvec = new Vector( tStop / dt )
	    Xgvec = new Vector( tStop / dt )
	    Xavec = new Vector( tStop / dt )
	    Xnvec = new Vector( tStop / dt )
	    for i = 0, 5 {
		tvec[i] = new Vector( tStop / dt )
		tvec[i].indgen( 0, dt )
	    }
	    gvec.record( &cell.syn[0].i )
	    avec.record( &cell.syn[1].i )
	    nvec.record( &cell.syn[2].i )
	    Xgvec.record( &cell.syn[3].i )
	    Xavec.record( &cell.syn[4].i )
	    Xnvec.record( &cell.syn[5].i )
	}
	vPlot.size( 150, 350, -70, -40 )
	run( pl, disp )
	cvec = gvec.c.add( avec ).add( nvec ).add( Xavec ).add( Xgvec ).add( Xnvec )
	tvec[1] = tvec[0].c.indvwhere( "<=", ns.start )
	low = tvec[1].x[ tvec[1].size()-1 ]
	/* Current is nA, time is ms, convert to density charge in uC/cm2
	* Cell area = 1e-6 cm2
	* 1nA * 1ms / 1AREA = 1pC/1e-6cm2 = 1uC/cm2 
	* ---> no conversion needed */
	charge = integral( tvec[0], cvec, low )
	gabaChrg = integral( tvec[0], gvec, low )
	ampaChrg = integral( tvec[0], avec, low )
	nmdaChrg = integral( tvec[0], nvec, low )
	XgabaChrg = integral( tvec[0], Xgvec, low )
	XampaChrg = integral( tvec[0], Xavec, low )
	XnmdaChrg = integral( tvec[0], Xnvec, low )
	
	ns = new Vector()
	return charge
    }
    
    /* Produce a PSP train. Records the injected charge.
    *
    * Arg 1, exin    : 0 = GABA synapse
    *                  1 = AMPA synapse
    *                  2 = NMDA synapse
    * Arg 2, w       : The synaptic weight in mS/cm2
    * Arg 3, interval: Time between PSPs in ms
    * Arg 4, nSpike  : The number of PSPs
    *(Arg 5, ext     : -1 = both local and external
    *                   0 = local
    *                   1 = external
    *(Arg 6, plot    : "plot" or "no plot")
    */
    proc PSPtrain() { local mode, interval, loc, exin, w, i
	exin = $1
	w = $2
	interval = $3
	nSpike = $4
	if( numarg() > 4 ) {
	    ext = $5
	} else {
	    ext = -1
	}
	if( numarg() > 5 && plotPoss ) {
	    pl = $s6
 	} else {
	    pl = "no plot"
	}
	delay = 0
	cell.disconnectCell()
	cell.soma ns = new MyNetStim(.5)
	ns.interval = interval
	ns.number = nSpike
	ns.start = 0
	ns.noise = 0
	print "interval ", interval, "number ", nSpike
	if( ext ) {
	    cell.pre_list[exin+3].append( new NetCon( ns, cell.syn[exin+3], 0, delay, 0 ) )
	    cell.setWsyn( exin, w, cvec, 0, 1 ) //w = 15
	}
	if( ext < 1 ) {
	    cell.pre_list[exin].append( new NetCon( ns, cell.syn[exin], 0, delay, 0 ) )
	    cell.setWsyn( exin, w )
	}
	cell.activate_syn( 1 )
	tStop = nSpike * interval + ns.start 
	if( cvode.active() ) {
	    tvec[0] = new Vector(10000)
	    tvec[1] = new Vector(10000)
	    Xavec = new Vector(10000)
	    avec = new Vector(10000)
	    cvode.record( &cell.syn[exin+2].i, Xavec, tvec[0] )	    
	    cvode.record( &cell.syn[exin].i, avec, tvec[1] )
	} else {
	    tvec[0] = new Vector( int( tStop / dt * 1.05 ) )
	    tvec[1] = new Vector( int( tStop / dt * 1.05 ) )
	    cvec = new Vector( int( tStop / dt * 1.05 ) )
	    avec = new Vector( int( tStop / dt * 1.05 ) )
	    Xavec = new Vector( int( tStop / dt * 1.05 ) )
	    avec.record( &cell.syn[1].i )
	    Xavec.record( &cell.syn[3].i )
	    tvec[0].indgen( dt )
	    tvec[1].indgen( dt )
	}
	if( strcmp( pl, "plot" ) == 0 ) {
	    vPlot.size( 0, tStop, -80, 20 )
	    vPlot2.size( 0, 10, 0, 0.001 )
	    spikevec = new Vector()
	    cell.pre_list[exin+3].object(0).record( spikevec )
	    run("plot")
	    for i = 0, spikevec.size()-1 {
		vPlot2.mark( spikevec.x[i], 0.0006, "o", 2 )
	    }
	    vPlot2.flush()
	    doEvents()
	} else {
	    run("no plot")
	}
	cvec = avec.c.add( Xavec )
	tvec[0] = tvec[0].c( 0, cvec.size()-1 )
	low = 0 /* Index for start of integral */
	/* Current is nA, time is ms, convert to density charge in uC/cm2
	* Cell area = 1e-6 cm2
	* 1nA * 1ms / 1AREA = 1pC/1e-6cm2 = 1uC/cm2 
	* ---> no conversion needed */
	curr = integral( tvec[0], cvec, low ) / ( (tvec[0].x[tvec[0].size()-1] - tvec[0].x[0] )/1000)
	/* 1uC/cm2 / (1ms/1000) = 1uC/cm2 / 1s = 1uA/cm2 */
	print "weight: ", w
    }
    
    /* VERSION NO 2 */
    /* Produce a PSP train
    *
    * Arg 1, exin    : 1 = AMPA synapse
    *                  2 = NMDA synapse
    * Arg 2, w       : The synaptic weight in mS/cm2
    * Arg 3, rate    : Rate of PSPs in Hz
    * Arg 4, tStop   : End of simulation
    *(Arg 5, plot    : "plot" or "no plot")
    */
    proc PSPtrain2() { local mode, interval, loc, exin, w, i
	w = $1
	rate = $2
	tStop = $3
	cell.disconnectCell()
	cell.poissonExternal(rate, w)
	if( numarg() > 3 && plotPoss ) {
	    pl = $s4
	} else {
	    pl = "no plot"
	}
	cell.activate_syn( 1 )
	//tStop = tStop + 150	
	if( strcmp( pl, "plot" ) == 0 ) {
	    vPlot.size( 0, tStop, -80, 20 )
	    vPlot2.size( 0, 100, 0, 0.001 )
	    spikevec = new Vector()	    
	    cell.pre_list[4+12*celltype].object(0).record( spikevec )
	    run("plot")
	    for i = 0, spikevec.size()-1 {
		vPlot2.mark( spikevec.x[i], 0.0006, "o", 2 )
	    }
	    vPlot2.flush()
	    doEvents()
	} else {
	    run("no plot")
	}
    }
    
    /* Produce several PSP trains. All arguments are vectors with 
    * parameters for the PSP trains, such that all parameters with a common 
    * vector index belong to the same PSP train. Parameters as below.
    *
    * Arg 1, start    : Start of netstim arrivals
    * Arg 2, del      : Delay before starting to measure activity
    * Arg 3, exinvec  : 1 if the channel is excitatory, 0 if it is inhibitory)
    * Arg 4, intervec : Time between PSPs in ms
    * Arg 5, nspikevec: The number of PSPs
    * Arg 6. shiftvec : The delay to the following PSP train
    * Arg 7. noisevec : The delay to the following PSP train
    *(Arg 8, wvec     : The weight of the synapse in mS/cm2)
    *(Arg 9, samevec  : If one wants several PSPs to come from the same source)
    * Return freq: Firing frequency in Hz
    */
    func PSPtrains() { local i, j, n1, n2, start, del, loc, freq
	freq = 0
	cell.disconnectCell()
	start = $1
	del = $2
	exinvec = $o3
	intervec = $o4
	nspikevec = $o5
	shiftvec = $o6
	noisevec = $o7
	if( numarg() > 7 ) {
	    wvec = $o8
	} else {
	    wvec = new Vector( intervec.size(), 0.1 )
	}
	if( numarg() > 8 ) {
	    samevec = $o9
	} else {
	    samevec = new Vector( intervec.size(), -1 )
	}
	objref nss[ intervec.size() ]
	for i = 0, intervec.size() - 1 {
	    loc = int( exinvec.x[i] / 6 ) 
	    if( i>0 && samevec.x[i] >= 0 ) {
		nss[ i ] = nss[ samevec.x[i] ]
	    } else {
		if( celltype == 1 ) {
		    cell.soma[ loc ] nss[ i ] = new MyNetStim(.5)
		} else {
		    cell.soma nss[ i ] = new MyNetStim(.5)
		}
	    } 
	    if( samevec.x[i] < 0 ) { /* Otherwise, this is already set */
		nss[ i ].interval = intervec.x[i]
		nss[ i ].number = nspikevec.x[i]
		if( i > 0 ) {
		    nss[ i ].start = start + shiftvec.x[i-1]
		} else {
		    nss[ i ].start = start
		}
		nss[ i ].noise = noisevec.x[i]
	    }
	    if( celltype == 1 ) {
		exin = ( exinvec.x[i] - 6 * loc )
		cell.connect_pre( nss[ i ], 0, exin, loc, 1 )
		cell.setWsyn( exin, wvec.x[i], loc )
	    } else {
		cell.connect_pre( nss[ i ], 0, exinvec.x[i], 1 )
		cell.setWsyn( exinvec.x[i], wvec.x[i] )
	    }
	}
	cell.activate_syn( 1 )
	//resetWvec()
	tvec = new Vector( intervec.size() )
	for i = 0, intervec.size() - 1 {
	    tvec.x[i] = intervec.x[i] * nspikevec.x[i]
	}
	tStart = start + del
	tStop = tvec.max() + tStart - del
	if( plotPoss ) {
	    vPlot.size( tStart, tStop, -80, 90 )
	    vPlot2.size( tStart, tStop, 0, 0.002 )
	    objref spikevec[ intervec.size() ]
	    for i = 0, intervec.size() - 1 {
		spikevec[i] = new Vector()	    
		cell.pre_list[ exinvec.x[i] ].object(0).record( spikevec[i] )
	    }
	    run("plot", 1, tStart )
	    for i = 0, intervec.size() - 1 {
		for j = 0, spikevec[i].size()-1 {
		    vPlot2.mark( spikevec[i].x[j], 0.0006*i, "o", 2 )
		}
	    }
	    vPlot2.flush()
	    doEvents()
	} else {
	    run("no plot", 1, tStart )
	}
	freq = nspk * 1000 / (tStop - tStart)
	return freq
    }
    
    /* Produces one spike 
    */
    proc oneSpike() { local amp, del, dur
	amp = 4
	del = 150
	dur = 6
	cell.setIapp( amp, del, dur )
	tStop = 500
	if( plotPoss ) {
	    vPlot.size( 100, 200, -70, 135 )
	    vPlot2.size( 155, 165, -15, 5 )
	    run("plot")
	} else {
	    run("no plot") 
	}
    }
    
    /* IClamp
    * Arg 1, amp: Amplitude in uA/cm2
    * Arg 2, dur: Duration in ms
    * Arg 3, str: no plotting if "no plot", otherwise plotting
    */
    proc iClamp() { local amp, dur, loc
	amp = $1
	del = 350
	dur = $2
	if( celltype == 1 ) {
	    loc = 2 /* Input only into distal dendrite */
	    cell.setIapp( amp, del, dur, loc )
	} else {
	    cell.setIapp( amp, del, dur )
	}
	tStop = del + dur + 100
	if( strcmp( $s3, "plot" ) == 0 && plotPoss ) {
	    vPlot.size( del, del+dur+30, -70, 70 )
	    vPlot2.size( del, del+250, -50, 25 )
	    //vPlot.size( del - 50, tStop, -80, 50 )
	    run($s3)
	} else {
	    run("no plot")
	}
    }
    
    proc resetWvec() { local i, j, k, m
	print "Resetting weight vector"
	for j = 0, cell.synCnt - 1 {
	    for k = 0, cell.pre_list[ j ].count() - 1 {
		for m = 1, cell.pre_list[ j ].object(k).wcnt() - 1 {
		    cell.pre_list[ j ].object(k).weight[m] = 0
		}
	    }
	}
    }
    
    /* A common function used by all other functions to run a simulation */
    proc run() { local i, disp
		
	/* Do simulation */
	if( numarg() > 0 && strcmp( pl1, "plot" ) == 0 ) {
	    pl = $s1
	} else {
	    pl = "no plot"
	}
	if( numarg() > 1 ) {
	    disp = $2
	} else {
	    disp = 1
	}
	if( numarg() > 2 ) {
	    tStart = $3
	} else {
	    tStart = 0
	}
	if( disp ) {
	    print ""
	    print "ECellTest.run"
	}
	t = 0
	nspk = 0
	steps = 0
	if( strcmp( pl, "plot" ) == 0 ) {
	    vPlot.erase()
	    vPlot2.erase()
	    reGraph()
	}
	finitialize()
	
	while( t < tStop ) {
	    if( cvode.active() ) {
		for i = 0, 9 {
	 	    fadvance()
		}
		steps = steps + 1
	    } else {
 		for i = 0, 9 {
	 	    fadvance()
		    steps = steps + 1
 		}
	    }
 	    if( cell.spike() == 1 && t > tStart ) {
	 	nspk = nspk + 1
	    }
	    
	    if( strcmp( pl, "plot" ) == 0 ) {
		vPlot.plot(t)
		vPlot2.plot(t)
	    }
	}
	if( strcmp( pl, "plot" ) == 0 ) {
	    vPlot.flush()
	    vPlot2.flush()
 	    doEvents()
	}
	if( disp ) {
	    print "finish time: ", t
	    print "ECellTest.run...  nspk: ", nspk
	    print "ECellTest.run... steps: ", steps
	}
	if( cvode.active() && printStat ) {
	    cvode.statistics()
	}
    }
    
    /* A general plot. Vectors xx and yy are plotted against eachother */
    proc plotXY() {
	siz = $o1.size()
	xx = $o1
	yy = $o2
	xmin = xx.min()
	ymin = yy.min()
	xmax = xx.max()
	ymax = yy.max()
	xyPlot.size( xmin, xmax, ymin, ymax ) 
	xyPlot.beginline()
	for i = 0, siz - 1 {
    	    xyPlot.line( xx.x[i], yy.x[i] )
	}
	xyPlot.flush()
	doEvents()
    }
    
endtemplate LabCell