demo_bombardment_passve_precise.oc

/*----------------------------------------------------------------------------

	CURRENT-CLAMP SIMULATIONS OF CORTICAL PYRAMIDAL CELLS
	=====================================================

	"precise" simulation: about 20,000 synapses simulated
	"passive": only passive properties

	Morphology
		- reconstructed Layer VI pyramidal cell from
		  Contreras, Destexhe and Steriade, 1997
		- correction for spines: 45% of dendritic membrane area
		- simple axon

	Passive properties
		- passive parameters adjusted to recordings in the absence
		  of synaptic activity (TTX + synaptic blockers)
		- passive parameters adjusted by simplex fitting to both
		  somatic and dendritic recordings
		  (dendritic recording: cell x210x4, Rin of 154 Meg after NBQX)
		=> Rin of 58.942 Meg in soma and 146 meg in dend1[12](0.179)
		
	Synaptic coverage:
		- AMPA and NMDA synapses in dendrites only; GABAa everywhere
		- exact synapse coverage for exc synapses (1.7um2)
		- exact synapse coverage for inh synapses (10um2)
		=> synapse densities consistent with morphological estimates
		   (DeFelipe & Farinas, 1992; Larkman 1991)

	Model adjusted to minis
		- quantal conductance compatible with patch-clamp (Sakmann)
		- uniform freq of release
		- parameters estimated from histograms
		=> gives minies with correct sigma and histograms

	Synaptic bombardment in passive model:
		- KCl: Erev_GABA = -55 mV; KAc: Erev_GABA = -75 mV
		- adjust freq to get 80% Rin change (Ketamine-Xylazine)
		- constrained by avg Vm under KCl (ECl=-55) and KAc (ECl=-75)

	Correlated bombardment:
		- correlated presynaptic random generator (corrGen8)

	Optimized algorithm:
		- multisynapse mechanisms in each segment
		=> tremendous acceleration of computation time


        Details of the models can be found in:

	Destexhe, A. and Pare D.  Impact of network activity on the integrative 
	properties of neocortical pyramidal neurons in vivo. J. Neurophysiol.
	81: 1531-1547, 1999.

	A PDF copy of this paper is available in http://cns.iaf.cnrs-gif.fr


	Alain Destexhe, destexhe@iaf.cnrs-gif.fr

----------------------------------------------------------------------------*/



//----------------------------------------------------------------------------
//  load and define general graphical procedures
//----------------------------------------------------------------------------

load_file("nrngui.hoc")

objectvar g[20]			// max 20 graphs
ngraph = 0

proc addgraph() { local ii	// define subroutine to add a new graph
				// addgraph("variable", minvalue, maxvalue)
	ngraph = ngraph+1
	ii = ngraph-1
	g[ii] = new Graph()
	g[ii].size(tstart,tstop,$2,$3)
	g[ii].xaxis()
	g[ii].yaxis()
	g[ii].addvar($s1,1,0)
	g[ii].save_name("graphList[0].")
	graphList[0].append(g[ii])
}

proc addshape() { local ii	// define subroutine to add a new shape
				// addshape()
	ngraph = ngraph+1
	ii = ngraph-1
	g[ii] = new PlotShape()
	g[ii].scale(-130,50)
}

nrnmainmenu()			// create main menu
nrncontrolmenu()		// create control menu




//----------------------------------------------------------------------------
//  transient time
//----------------------------------------------------------------------------

CURRINJ = 0		// amount of injected current - serves as flag

if(CURRINJ == 0) {
  trans = 150	// transient to reach steady state
} else {
  trans = 300	// transient to skip injected current
}

v_init = -65		// initial condition

print " "
print ">> Transient time of ",trans," ms"
print " "


DEBUG=0


//----------------------------------------------------------------------------
//  create multi-compartment geometry
//----------------------------------------------------------------------------

print " "
print ">> Reading geometry of neuron..."
print " "

xopen("layer6.geo")		// Layer VI pyramidal cell

corrD = 1.449		// dendritic correction for spines (44% of membrane)



//----------------------------------------------------------------------------
//  add a simple axon
//----------------------------------------------------------------------------

xopen("add_just_axon.oc")		// add simplified axon




//----------------------------------------------------------------------------
//  Passive currents
//----------------------------------------------------------------------------

// Best fit for TTX-bicuculline with Layer 6 cell, soma
// fixed: rev=-65, cm=1, Ra=250, corrD=1.449
// fit: g_pas=4.52e-5  (Error=5.0221802)

leak_cond = 4.52e-5
leak_rev = -65
leak_rev = -70			// adjusted to cell x210x4
leak_rev = -80			// fr3
capacit = 1
axial_res = 250

forall { 			// insert passive current everywhere
	insert pas
	g_pas = leak_cond
	e_pas = leak_rev
	cm = capacit
	Ra = axial_res
	L = L
}

forsec "axon" {  		// exceptions along the axon
	cm = 0.04
	g_pas = 0.02
}

forsec "dend" { 		// correction for dendrites
	g_pas = g_pas * corrD
	cm = cm * corrD
}




//----------------------------------------------------------------------------
//  localize synapses
//----------------------------------------------------------------------------

// Nov 27, 1997: recalculated densities to make them compatible with the 
// proportion of synapses found in pyramidal cells

cutoff = 40		// cutoff distance (um) where spines begin
ex_dend_unit =  1.7  // 100 // unit membrane area for excitatory synapses
in_dend_unit =  10   // 100 // unit membrane area for inh synapses in dendrites
in_soma_unit =  2.5  // 25  // unit membrane area for inh synapses in soma
in_iseg_unit =  1.7  // 17  // unit membrane area for inh synapses in init seg

//  With 100,100,25,17 um2 (exc dend, inh dend, inh soma, inh iseg), one
//  excitatory synapse represents 55-65 real synapses and one inhibitory
//  synapse represents 8.8-10.4 real synapses... (ratio of 6.25)
//  (according to high spine density; and 7% GABAergic in soma)


xopen("localize_synapses_corrgen_mul.oc")   // procedures and initializations

SEED = 1				  // flag for seed
if(SEED) set_seed(0.1,0.2,0.3,0.4)	  // set seed for random numbers


EXC = 1		// flag variable to insert excitatory synapses
NMDA = 0	// flag variable for NMDA
INH = 1		// flag variable to insert inhibitory synapses

if(INH) {
  insert_GABA_prox()	// insert GABAa synapses in soma, prox dend & axon
  insert_GABA_dend()	// insert GABAa synapses in dendrites
}

if(EXC) { 
  insert_AMPA_dend()	// insert AMPA synapses in dendrites
  if(NMDA) {
     insert_NMDA_dend()	// insert NMDA synapses in dendrites
  }
}



//
//  Presynaptic parameters
//
pre_freq_I = 5.5 	// inh presynaptic frequency
pre_freq_E = 1.0 	// exc presynaptic frequency 
			// (if inh is 0.1, exc should be 0.625)
pre_dur = 1e6		// duration of presynaptic firing

corr_E = 0.7	  	// exc correlation
corr_I = 0.7	  	// inh correlation

set_generators()



//
// KINETICS
//

Erev_multiGABAa = -55	// chloride (from Denis)
Erev_multiGABAa = -75	// K-Ac

//Cdur_multiGABAa = 0.3
//Alpha_multiGABAa = 20
//Beta_multiGABAa = 0.05	// from SimFit to Denis recordings
//Beta_multiGABAa = 0.18	// SimFit to hippocampal GABAa

Cdur_multiGABAa = 1	// idem Meth Neuronal Modeling
Cmax_multiGABAa = 1	// idem Meth Neuronal Modeling
Alpha_multiGABAa = 5	// idem Meth Neuronal Modeling
Beta_multiGABAa = 0.1	// fr3

//Cdur_multiAMPA = 0.3
//Alpha_multiAMPA = 5
//Alpha_multiAMPA = 20	// better (higher amplitude)
//Beta_multiAMPA = 0.243	// from SimFit to Denis recordings

Cdur_multiAMPA = 1		// idem Meth Neuronal Modeling
Cmax_multiAMPA = 1		// idem Meth Neuronal Modeling
Alpha_multiAMPA = 1.1	// idem Meth Neuronal Modeling
Beta_multiAMPA = 0.67	// fast AMPA to get a decay of 1.5 ms (Markram)



//
// QUANTAL CONDUCTANCES
//

g_AMPA = 0.001200	// quantal AMPA conductance (Denis is 0.000260)
g_GABA = 0.000600	// quantal GABA conductance (consistent with in vitro)
 
// By comparison, Sakmann is 200-400 nS for GABA, AMPA is 0.35-1 nS (McBain
// and Dingledine, 1992; Burgard and Hablitz, 1993)

if(EXC) {
  if(NMDA) {
     g_NMDA = 4 * g_AMPA
  } else {
     g_NMDA = 0
  }
} else {
  g_AMPA = 0
  g_NMDA = 0
}

if(INH) { 
  		// do nothing
} else {
  g_GABA = 0
}






proc stim_uniform() {
  set_generators()
  if(EXC) {
     set_AMPA_dend(g_AMPA*corrD)	// dendritic AMPA conductances
     if(NMDA) {
	set_NMDA_dend(g_NMDA*corrD)	// dendritic NMDA conductances
     } 
  }
  if(INH) {
     set_GABA_prox(g_GABA)		// perisomatic GABA conductances
     set_GABA_dend(g_GABA*corrD)	// dendritic GABA conductances
  }
  printf("\nSetting generators and synaptic conductances:\n")
  printf(" Exc f = %g Hz\n Inh f = %g Hz\n",pre_freq_E,pre_freq_I)
  printf(" gAMPA = %g uS\n gNMDA = %g uS\n gGABA = %g uS\n", \
	g_AMPA,g_NMDA,g_GABA)
}

stim_uniform()




//----------------------------------------------------------------------------
//  insert electrode in dendrite or soma
//----------------------------------------------------------------------------

xopen("Electrode.oc")		// template for electrode

access soma
//access dend1[12]

objectvar El		// create electrode

El = new Electrode(0.5)

soma El.stim.loc(0.5)		// locate in soma
//dend1[12] El.stim.loc(0.179)	// locate in dendrite
El.stim.del = 0
El.stim.dur = 1e6
El.stim.amp = 0



objectvar dc		// create DC-current

dc = new Electrode(0.5)

soma dc.stim.loc(0.5)		// locate in soma
//dend1[12] dc.stim.loc(0.179)	// locate in dendrite
dc.stim.del = 0
dc.stim.dur = 1e6
dc.stim.amp = 0





//----------------------------------------------------------------------------
//  setup simulation parameters
//----------------------------------------------------------------------------

Dt = 0.1
npoints = 10000	// 600000

objectvar SIMsoma,SIMdend		// create vectors of simulation points
SIMsoma = new Vector(npoints+500)
SIMdend = new Vector(npoints+500)

dt = 0.1			// must be submultiple of Dt
tstart = 0
tstop = npoints * Dt
runStopAt = tstop
steps_per_ms = 1/Dt
celsius = 36

statpts = npoints+1-trans/Dt	// nb of points to analyze

objectvar Vsoma, Vdend		// create vectors for histogram analysis
Vsoma = new Vector(statpts)
Vdend = new Vector(statpts)





//----------------------------------------------------------------------------
//  Define histogram procedures
//----------------------------------------------------------------------------

nbins = 100			// nb of points in histogram
vmin = -80			// min value of Vm
vmax = 0			// max value of Vm
hmax = 20000			// max value of histogram
binsize = (vmax-vmin)/nbins	// size of bin

objectvar Hsoma,Hdend		// create vectors for histograms
Hsoma = new Vector(nbins)
Hdend = new Vector(nbins)
objectvar HX
HX = new Vector(nbins)		// Vector for histogram's absissa
x = vmin
for i=0, nbins-1 {
	HX.set(i,x)
	x = x + binsize
}

hgr = ngraph
g[hgr] = new Graph()		// graph for histogram
g[hgr].size(vmin,vmax,0,hmax)
g[hgr].xaxis()
g[hgr].yaxis()
g[hgr].save_name("graphList[0].")
graphList[0].append(g[hgr])
ngraph = ngraph + 1

proc init() {				// initialization procedure
	finitialize(v_init)
	fcurrent()
	index = 0			// add definition of an index
}


proc step() {local i			// advance-one-step (Dt) procedure
	Plot()
	SIMsoma.set(index,soma.v(0.5))		// memorize data
	SIMdend.set(index,dend1[12].v(0.179))	// memorize data
	index = index + 1
	for i=1,nstep_steprun {
		advance()
	}
}

//
//  calculate sigma from histogram (skipping spikes)
//
proc calc_sigma() { local sum,avg,sig
   x = vmin
   for i=0, nbins-1 {
	if(x <= $1) {
		y = Hsoma.get(i)
		sum = sum + y
		avg = avg + y * x
		sig = sig + y * x*x
	}		
	x = x + binsize
   }
   avg = avg / sum
   sig = sqrt(sig/sum - avg*avg)
   printf("\n=> Values computed by cutting spikes: avg=%g, sigma=%g\n\n",avg,sig)
}

niter = 1
Rin = 0

proc run_histo() {
   for i=0, niter-1 {
	if(SEED) set_seed(0.1,0.2,0.3,0.4)		// set seed
	run()						// run simulation

	Vsoma.copy(SIMsoma,trans/Dt,npoints-1)		// truncate data
	Hsoma = Vsoma.histogram(vmin,vmax,binsize)	// make histogram
	Hsoma.plot(g[hgr],HX)				// draw histogram
	Avg = SIMsoma.mean(trans/Dt,npoints-1)		// calc statistics
	Std = SIMsoma.stdev(trans/Dt,npoints-1)
	if(CURRINJ != 0) {
	 Rin=-(SIMsoma.mean(320/Dt,400/Dt)-SIMsoma.mean(120/Dt,200/Dt))/CURRINJ
	}
	printf("\nSoma:\tRin=%g\tAvg=%g\tStd=%g\n",Rin,Avg,Std)
	calc_sigma(-40)

	Vdend.copy(SIMdend,trans/Dt,npoints-1)		// truncate data
	Hdend = Vdend.histogram(vmin,vmax,binsize)	// make histogram
	Hdend.plot(g[hgr],HX)				// draw histogram
	Avg = SIMdend.mean(trans/Dt,npoints-1)		// calc statistics
	Std = SIMdend.stdev(trans/Dt,npoints-1)
	if(CURRINJ != 0) {
	 Rin=-(SIMsoma.mean(320/Dt,400/Dt)-SIMsoma.mean(120/Dt,200/Dt))/CURRINJ
	}
	printf("dend:\tRin=%g\tAvg=%g\tStd=%g\n",Rin,Avg,Std)
   }
}



proc make_SBpanel() {			// make panel
	xpanel("Syn Bombardment")
	xpvalue("g_AMPA",&g_AMPA)
	xpvalue("g_NMDA",&g_NMDA)
	xpvalue("g_GABA",&g_GABA)
	xpvalue("Exc freq",&pre_freq_E)
	xpvalue("Inh freq",&pre_freq_I)
	xpvalue("Exc correlation",&corr_E)
	xpvalue("Inh correlation",&corr_I)
	xpvalue("Cl reversal",&Erev_multiGABAa)
	xpvalue("AMPA decay",&Beta_multiAMPA)
	xpvalue("GABA decay",&Beta_multiGABAa)
	xbutton("Apply","stim_uniform()")
	xbutton("Set seed","set_seed(0.1,0.2,0.3,0.4)")
	xpvalue("Nb iterations",&niter)
	xbutton("Run + calc histogram","run_histo()")
	xpanel()
}

make_SBpanel()




//----------------------------------------------------------------------------
//  add graphs
//----------------------------------------------------------------------------

addgraph("soma.v(0.5)",vmin,vmax)		// soma
addgraph("dend1[12].v(0.179)",vmin,vmax)