/* De Schutter model conductances in De Schutter-Rapp cell      */
/* conversion from GENESIS by Jenny Davie, Arnd Roth,
   Volker Steuber, Erik De Schutter & Michael Hausser 28.8.2004 */

/********************************************/
/************ initial parameters ************/
/********************************************/

t = 0                  /* simulation starts at t = 0 ms */
dt = 0.02              /* time step, ms */
Vrest = -68            /* resting potential, mV */
celsius = 37           /* invoke Q10 for conductances and kinetics */
tstop = 1300           /* simulation stops at t = 1300 ms */

strdef outputfile, prefix
prefix = "DS2M0PurkVar2"

/********************************************/
/*************** make geometry **************/
/********************************************/

xopen("Purk2M0.nrn")   /* load the morphology file */

/*************** passive membrane properties  ***********/

membranecap = 1.64     // specific membrane capacitance in uF cm^-2
membraneresist = 30000 // specific membrane resistance in ohm cm^2 for dendrite
somamembraneresist = 10000
axialresist = 250      // axial resistivity in ohm cm

forall {insert Leak el_Leak=-80 Ra=axialresist}

forsec "dend" {gl_Leak=1.0/membraneresist Ra=axialresist cm=1.0*membranecap}

forsec md {gl_Leak=1.0/membraneresist Ra=axialresist cm=1.0*membranecap}

forsec "soma" {gl_Leak=1/somamembraneresist Ra=axialresist cm=1.0*membranecap}

spine_dens = 13
spine_area = 1.33

func spinecomp() { local a
	maxdiam = 0
	for (x) if (diam(x) > maxdiam) maxdiam = diam(x)
	if (maxdiam <= 3.17) { // spine correction only for thin dendrites
		a = 0
		for(x) a=a+area(x)

		F = (L*spine_area*spine_dens + a)/a
	} else {
		F = 1 // no spine correction for larger dendrites
	}
	return F
}

proc add_spines() {
		spinecomp()
		gl_Leak = spinecomp()/membraneresist
		cm = spinecomp()*membranecap
}

forsec "dend" add_spines()


/*************** insert active conductances  ***********/

shell_depth = 0.2 // um

soma {
insert NaF gnabar_NaF = 7.5
insert NaP gnabar_NaP= 0.001 // default is correct value for soma
//insert CaP gcabar_CaP = 0.0 // no P type Ca Ch in soma
insert CaT gcabar_CaT = 0.0005 // default is correct value for soma and dend
insert Kh1 gkbar_Kh1 = 0.0003 // default is correct value for soma, Ih uses only K ions here
insert Kh2 gkbar_Kh2 = 0.0003 // default is correct value for soma, Ih uses only K ions here
insert Kdr gkbar_Kdr = 0.6 // default is correct value for soma
insert KMnew2 gkbar_KMnew2 = 0.00004 // default is correct value for soma
insert KA gkbar_KA = 0.015 // default is correct value for soma
//insert KC gkbar_KC = 0.0 // no BK Ch in soma
//insert K2 gkbar_K2 = 0.0 // no KCa in soma
insert cad
equivalentDepth = shell_depth - shell_depth^2/diam
// GENESIS calculates shell area as (outer circle area - inner circle area)
// our calcium diffusion had used depth * outer circumference
depth_cad = equivalentDepth

cai = 4e-5 cao = 2.4
eca = 12.5*log(cao/cai)
ion_style("ca_ion", 1, 1, 0, 0, 0) // was ("ca_ion", 1, 1, 0, 1, 0)
ena = 45 // mV Na+ reversal potential
ek = -85 // mV K+ reversal potential
}

forsec "dend" {

/* gmax as in DeSchutter 1994 PM9 'rest of dendrite' */

//insert NaF gnabar_NaF = 0.0
//insert NaP gnabar_NaP = 0.0
insert CaP gcabar_CaP = 0.0045
insert CaT gcabar_CaT = 0.0005
//insert Kh1 gkbar_Kh1 = 0.0
//insert Kh2 gkbar_Kh2 = 0.0
//insert Kdr gkbar_Kdr = 0.0
insert KMnew2 gkbar_KMnew2 = 0.000013
//insert KA gkbar_KA = 0.0
insert KC gkbar_KC = 0.08
insert K2 gkbar_K2 = 0.00039
insert cad
equivalentDepth = shell_depth - shell_depth^2/diam
// GENESIS calculates shell area as (outer circle are - inner circle area)
// our calcium diffusion had used depth * outer circumference
depth_cad = equivalentDepth

cai = 4e-5 cao = 2.4
ion_style("ca_ion", 1, 1, 1, 1, 0) // was ("ca_ion", 1, 1, 0, 1, 0)
//ena = 45 // mV Na+ reversal potential
ek = -85 // mV K+ reversal potential
}

/* add (the same) conductences to main dend, the spine necks and spine heads*/

forsec md {
/* gmax as in DeSchutter 1994 PM9 'Main Dendrite'*/
//insert NaF gnabar_NaF = 0.0
//insert NaP gnabar_NaP = 0.0
insert CaP gcabar_CaP = 0.0045
insert CaT gcabar_CaT = 0.0005
//insert Kh1 gkbar_Kh1 = 0.0
//insert Kh2 gkbar_Kh2 = 0.0
insert Kdr gkbar_Kdr = 0.06
insert KMnew2 gkbar_KMnew2 = 0.000010
insert KA gkbar_KA = 0.002
insert KC gkbar_KC = 0.08
insert K2 gkbar_K2 = 0.00039
insert cad
equivalentDepth = shell_depth - shell_depth^2/diam
// GENESIS calculates shell area as (outer circle are - inner circle area)
// our calcium diffusion had used depth * outer circumference
depth_cad = equivalentDepth

cai = 4e-5 cao = 2.4
ion_style("ca_ion", 1, 1, 1, 1, 0) // was ("ca_ion", 1, 1, 0, 1, 0)
//ena = 45 // mV Na+ reversal potential
ek = -85 // mV K+ reversal potential
}

access soma
/*************** now insert current pulses and synapses *********/

	objectvar currentpulse

	currentpulse = new IClamp(0.5)       // in middle of soma
	currentpulse.del = 100 
	currentpulse.dur = 1000
	currentpulse.amp = 3.0


/********************************************/
/***** initialization of simulation run *****/
/********************************************/

proc run() {
print currentpulse.amp

	// initialization of simulation run
	t = 0                   // simulation starts at t = 0 ms
	count = 0
	finitialize(Vrest)
	fcurrent()

	while (t < tstop) {
		fadvance()
		if((count/5) == int(count/5)) {		
		fprint("%g %g %g %g %g\n", t, soma.v(0.5), dend[4].v(0.5), dend[19].v(0.5), dend[1511].v(0.5))
		}
	count += 1
	}

	wopen()
}

/******* generate an IV *********/

proc metarun() {local myvar

	for(myvar = 1; myvar < 2; myvar = myvar + 1) {

		scalar=1.5*myvar
		currentpulse.amp = scalar
		
		if (scalar<=0) {
		sprint(outputfile, "%sm%g", prefix, -scalar)
		}
		if (scalar>0) {
		sprint(outputfile, "%sp%g", prefix, scalar)
		}
		wopen(outputfile)
		run()
	}
}

metarun()

quit()