synresp.hoc

// For compatibility with MS Notepad, set tabs to every 8 char.

// Set major run control parameters (0=False=No, 1=True=Yes)
isCC=0				// set CC mode (otherwise VC)
isInteractive=1 		// run interactive vs in batch
runStim=0			// run the synaptic stimulation simulation

// Load the Neuron machinery
if (isInteractive) {
	load_file("nrngui.hoc")
} else {
	load_file("stdrun.hoc")
}
cvode.active(1) 		// use an adaptive method


// The following are set before this hoc file is invoked.
// Values shown are examples only.

// Define an orientation vector for locating layers.
// These component must be unit length and have
// point in the apical direction of the reoriented
// y coordinate for the cell.

// orientX=0			// No real orientation
// orientY=1			// so the unit vector is
// orientZ=0			// just the y axis direction.

// Locate the (reoriented) boundaries between layers
// PPy3d=350			// Minimum y-coord value for PP (PP-SR boundary)
// SRy3d=100			// Minimum Y-coord value for SR (SR-SL boundary)
// SOy3d=0			// Maximum y-coord value for SO


// Define input and output files (examples only)
// strdef cell			// hoc file with cell geometry
// strdef outPath		// Synaptic responses file name
// strdef avgPath		// Overall average response file name
// cell="bar-cell1zr.CNG.hoc"	// Cell geometry file
// outPath="out-temp.txt"	// Synaptic responses output file

// Provide the path for the common axon segment definition
strdef axonPath
axonPath="axon-common.hoc"

// Set a default path for saved response data vectors
strdef savePath
savePath="savedsynresp.csv"


// Simulation parameters -----------------------------

addChannels=0			// Add active channels to the cell
useHemond2008=0			// Use conductances from the Hemond et al., 2008 model.
useSpine=0			// Use a spine for synapses in SR and SO
maxSegLen=10			// Maximum segment length in microns

tError=1000 			// Maximum time after init before error declared (see runIt)
tInit=2000			// Initialization time to reach equilibrium
tCont=50			// Continue time (must be greater than time to peak)
tDelta=0.1			// Time step size for response simulation and averaging
randSeed=19720107		// Random number generator seed value
maxSecnum=1000			// maximum section number (see findLayer and onPPtrunk)

useAMPAR=1			// Run for AMPAR synapses (ie AMPAR not blocked)
useNMDAR=0			// Run for NMDAR synapses (ie NMDAR not blocked)

fastAMPAR=0			// use RC AMPAR synapse settings with same kinetics as PP
randomizeSyn=0			// make random variations to synaptic properties

maxDendDiam=5			// Maximum dendrite diameter (otherwise soma)

Vrest = -61 			// Starting resting potential everywhere
// Vrest = -61 - 13.6 		// Starting resting potential using LJP from Hemond et al.
// Vrest = -71 			// Starting resting potential shifted by 10 mV


celsius = 32.0			// Temperature (needed for active channels)
synErev = 0 			// Synaptic reversal potential (AMPAR and NMDAR)
lowExtMg = 50e-3		// External [Mg++] (mM) for low-Mg experiments
highExtMg = 1			// External [Mg++] (mM) for normal Mg experiments

VCvm = -80			// Voltage clamp potential at soma
CCvm = -60			// Current clamp holding potential - as per model paper predictions
// CCvm = -80			// Current clamp holding potential - more or less per experiments

// PP synapse parameters for AMPAR responses.
// Generally these are similar to values found in MF terminals
// and SC synapses in CA1 measured near the synapse. 
ppAMPARSynTau1=0.4		// PP synaptic rising time constant
ppAMPARSynTau2=4.1		// PP synaptic falling time constant
ppAMPARSynWeight=0.9		// PP synapse gmax in nS
ppAMPARSynTsd=0.0		// PP synapse tau2 log-normal sd param
ppAMPARSynWsd=0.2		// PP synapse weight log-normal sd param

// RC synapse parameters for AMPAR responses.
// These values are similar to values from Williams & Johnston 1991.
// The slower kinetics may result from combined AMPA/KA synapses,
// from reduced neurotransmitter concentration in the cleft,
// or even from multi-quantal release, all of which are hypothetical.
rcAMPARSynTau1=3.3		// RC synaptic rising time constant
rcAMPARSynTau2=3.3		// RC synaptic falling time constant
rcAMPARSynWeight=0.5		// RC synapse gmax in nS
rcAMPARSynTsd=0.0		// RC synapse tau2 log-normal sd param
rcAMPARSynWsd=0.4		// RC synapse weight log-normal sd param

if (fastAMPAR) {
	// Set RC AMPAR time constants to match fastest experimental pVC response
	// rcAMPARSynTau1=1.3
	// rcAMPARSynTau2=5.3

	// Set RC AMPAR to match PP AMPAR or other faster AMPAR model
	rcAMPARSynTau1=0.4
	rcAMPARSynTau2=4.1
}


// PP and RC synapse parameters for NMDAR responses.
// This is pure data fitting. The values are somewhat
// in line with Vicini et al. 1998 for HEK transfection,
// but there should be a larger correction for temperature.
// See also Dalby and Mody 2003 for NR2A in DG (similar).
// See Chen et al., 2001 Mol Pharamcol 59: 212-219 for one
// of the few (or only) experimental studies to consider 
// both NMDAR subtypes and temperature effects together.
// Mean NMDAR experimental values from Chen et al. are 
// 2.4 ms and 35 ms at 37 deg-C in outside-out patches. 
// For the parameters below, 10%-90% rise time is 4.7 ms while
// 90%-10% decay time is 37.9 ms. There is no good reason to
// assume the same gmax for RC and PP since the synapses
// have morphological differences.

ppNMDARSynTau1=5		// PP synaptic rising time constant
ppNMDARSynTau2=16		// PP synaptic falling time constant
ppNMDARSynWeight=0.18		// PP synapse gmax in nS
ppNMDARSynTsd=0.0		// PP synapse tau2 log-normal sd param
ppNMDARSynWsd=0.0		// PP synapse weight log-normal sd param
 
rcNMDARSynTau1=5 		// RC synaptic rising time constant
rcNMDARSynTau2=16		// RC synaptic falling time constant
rcNMDARSynWeight=0.16		// RC synapse gmax in nS
rcNMDARSynTsd=0.0		// RC synapse tau2 log-normal sd param
rcNMDARSynWsd=0.0		// RC synapse weight log-normal sd param 

// Mean passive values estimated in Hemond et al. 2009 (Ih paper)
// RaAll is not a sensitive fit in that context but a
// more typical value is needed to fit synaptic responses.
// For justification of different spine factors in PP and RC
// see Megias et al. 2001 for results in CA1 as well as 
// Matsuda et al. 2004 plus Ishizuka et al. 1995 for CA3.
// Axon properties are from Hemond et al. 2008

Cm=0.72 			// Capacitivity in microF/cm^2 (1.44/2)
Rm=62996			// Membrane resistivity in ohm-cm^2 (31498*2)

RaAll=140			// Axial resistivity in ohm-cm
//RaAll=198			// Axial resistivity in ohm-cm

RaAxon=50			// Axon-only axial resistivity
RaSpine=140			// Spine-only axial resistivity

PPSpineAdj=1.0			// Cm,Rm adjustment for PP spines
RCSpineAdj=2.0			// Cm,Rm adjustment for RC spines

// Parameters for hypothetical K+ channels in spine
ekspine = -90			// K+ reversal for channels in spine (in mV)
// gkspine = 120e-3		// K+ density actived in spine (in S/cm^2)
// gkspine = 55e-3		// K+ density actived in spine (in S/cm^2)
gkspine =0			// K+ density actived in spine (in S/cm^2)

// Provide channel conductances. These are in S/cm^2.
// In the model from Hemond et al., 2008 fig 9, there is a shift
// of voltage sensitivity by 24 mV in most channels.
// These conductance settings go with figure 9d (non-delay case)

gna=22e-3			// Sodium channel conductivity (all but axon)
gnaAxon=110e-3			// Sodium channel conductivity for axon
gkdr=10e-3			// K-dr (delayed rectifier) channel conductivity
gkap=20e-3			// K-A channel conductivity
gc=0.01e-3			// Ca-L,N,T channel conductivity in soma (all the same)
ghd=0.00e-3			// Ih channel conductivity 
gKc=0.05e-3			// K-C channel conductivity 
gahp=0.0e-3			// K-ahp channel conductivity 
gkm=0e-3			// K-M channel conductivity (adapting - soma only in this model)
gkd=0e-3			// Kd channel conductivity (delayed onset - soma only in this model)

sh=24				// Global voltage shift in mV for channels (Hemond et al. only)

gkaSlope=5.2/350		// Distance dependency for K-A (not used for Hemond et al.)
gkaSlopeMaxDist=9999		// Maximum distance for K-A dependency (not used for Hemond et al.)			


// Misc working variables (global) ----------------------------------
objref ns			// Common network stimulation
objref amparNC, nmdarNC		// Network connections
objref amparSyn, nmdarSyn	// Synapses
objref hold, seClamp		// Clamps
objref tvec, ivec, vvec		// Recording vectors for time, current, and voltage 
objref dvec			// Recording vector for dendrite Vm
objref iampar,inmdar		// Recording vector for AMPAR and NMDAR currents
objref outfile			// Output file
objref impedance		// For figuring out CC input impedance
objref saveFile			// File to save individual traces to (see saveSynResp below)
objref rand			// Random number generator 
objref strFun			// StringsFunction object
objref adOnPPpath		// flags lookup vector for apical_dendrites
objref secRef			// temporary section reference 

objref nil			// NULL object (so we can test for NULL, WTG hoc)

isError=0			// global error indicator		

strdef curSecname		// Name part of the current section (see parseSecname)
curSecnum=0			// Number part of the current section (see parseSecname)

nSomaSec=0			// Number of sections making up the soma
somaMidSec=0			// Index to section at the middle of the soma (roughly speaking)
strdef layer			// layer of current segment as string
layerId=0			// layer as number: 0=SOMA, 1=AXON, 2=SO, 3=SL, 4=SR, 5=LM, -1=ERROR

ory3d=0 			// reoriented y3d value for current segment
synLoc=0			// location of synapse in current section
segLen=0			// length of the current segment

ivec = new Vector()		// record of soma current injection
vvec = new Vector() 		// record of soma voltage
tvec = new Vector() 		// record of time
dvec = new Vector()		// record of voltage at the synapse
iampar = new Vector()		// record AMPAR current
inmdar = new Vector()		// record NMDAR current

adOnPPpath = new Vector(maxSecnum) // vector for onPPpath lookup by section number
strFun = new StringFunctions()	// functional instance for string functions


// We will accumulate some stats by layer etc.
// Trunk segments are those in SR that have a PP connection in their subtree.
totSomaLen=0			// Total length of all soma segments
totSomaArea=0			// Total area of all soma segments

totAxonLen=0			// Total length of all axon segments
totAxonArea=0			// Total area of all axon segments

totLMlen=0			// Total length of all segments in LM
totLMarea=0			// Total area of all segments in LM

totSRlen=0			// Total length of all segments in SR
totSRarea=0			// Total area of all segments in SR
totSRTlen=0			// Total length of trunk segments in SR 
totSRTarea=0			// Total area of all trunk segments in SR

totSOlen=0			// Total length of all segments in SO
totSOarea=0			// Total area of all segments in SO		

// Create spine sections
create spineHead // area will be roughly 1 micron^2
create spineNeck // resistance will be roughly 500 mega-ohm
spineHead {L=0.5 diam=0.6 nseg=1}
spineNeck {L=1 diam=0.06 nseg=1}
spineNeck connect spineHead(0),1


// Load the cell topology, otherwise
// Neuron gets confused about sections.
// A common axon segment is added below.
xopen(cell)   

// Remove any type 10 (marker) segments
if (section_exists("user10")) {
	print "Removing user 10 segments"
	forsec "user10" {
		disconnect()
		delete_section()
	}
}

// Remove any axon segments
if (section_exists("axon")) {
	print "Removing axon segments"
	forsec "axon" {
		disconnect()
		delete_section()
	}
}

// Now load a replacement axon segment
load_file(axonPath)


// ========================================================
// Subroutines (defined before they are referenced)
// ========================================================

// Utility routines for min/max of two values.
// Note that the documentation for proc even provides
// the code for this but there is no such min or max
// function outside of vectors.
func minVal() {if ($1<$2) { return $1 } else { return $2} }
func maxVal() {if ($1>$2) { return $1 } else { return $2} }


// Get the name and number of the current section. The
// section name is assumed to be in the form: name[num]
// Results are in globals curSecname and curSecnum.
proc parseSecname() {local n

	strdef stemp
	stemp=secname()
	n = strFun.substr(stemp,"[")
	curSecname=stemp
	curSecnum=-1
	if (n>=0) {
		strFun.left(curSecname,n)
		strFun.right(stemp,n)
		sscanf(stemp,"[%d",&curSecnum)
	}
}
	
	
// Find the layer associated with the current segment.
// xseg is the relative location within the section (0 to 1).
// ory3d is set based on segment location in the reoriented
// cell morphology using orientX, orientY, and orientZ.
// layer is set to one of "SOMA", "AXON", "SO", "SL", "SR", 
// "LM", or "ERROR" with layerId set accordingly (0-5,-1). 
// layerId is returned because you can set string values
// in hoc but cannot test their values except via strcmp.
// onPPpath is set to 1 if the current section is an
// apical_dendrite with a PP segment in its subtree.
// Otherwise onPPpath is 0. curSecname and curSecnum are
// set as a side-effect of calling parseSecname.

proc findLayer() {local xseg,xarc,ipt3d,al0,al1,sx3d,sy3d,sz3d

	xseg = $1
	layer = "ERROR" // no matching layer found (so far)
	layerId=-1
	onPPpath=0

	// Determine if the current section is an
	// apical_dendrite with a PP segment in its
	// subtree. adTrunk vector is assumed to
	// already be set with the appropriate values
	parseSecname()
	if (strcmp(curSecname,"apical_dendrite")==0) {
		onPPpath=adOnPPpath.x(curSecnum)
	}

	// If the current section is a spineHead or spineNeck,
	// assign an (so far unused) layer and id and stop
	if (strcmp(curSecname,"spineNeck")==0 || strcmp(curSecname,"spineHead")==0) {
		layer="SPINE"
		layerId = 999
		return
	}

	// Check for maximum allowed diameter.
	// Larger dendrites are most likely to be
	// really part of the soma regardless of type.
	if (diam(xseg)>maxDendDiam || issection("soma.*")) {
		layer = "SOMA"
		layerId = 0
		return
	}

	// Also look for any axon segments just in case
	if (issection("axon.*")) {
		layer = "AXON"
		layerId = 1
		return
	}

	// Locate 3D points for the current segment
	// to get a Y-coordinate value.
	xarc=xseg*L
	ipt3d=1
	while (ipt3d<n3d() && arc3d(ipt3d)<xarc) { ipt3d += 1 }
	if (ipt3d<n3d()) {

		// Estimate the segment center location
		// from a linear interpolation of points
		// before and after it in arc length.
		al0=xarc-arc3d(ipt3d-1)
		al1=arc3d(ipt3d)-xarc
		sx3d=(al0*x3d(ipt3d)+al1*x3d(ipt3d-1))/(al0+al1)
		sy3d=(al0*y3d(ipt3d)+al1*y3d(ipt3d-1))/(al0+al1)
		sz3d=(al0*z3d(ipt3d)+al1*z3d(ipt3d-1))/(al0+al1)

		// Get the distance along the reoriented y-coord
		// orientX,Y,Z form a unit vector in the proper
		// direction with the apical direction up (>0).
		ory3d=sx3d*orientX+sy3d*orientY+sz3d*orientZ

		// Classify layer based on oriented y-coord
		if (ory3d>=PPy3d) {
			layer="LM"
			layerId=5
		} else if (ory3d>=SRy3d) {
			layer = "SR"
			layerId=4
		} else if (ory3d>SOy3d) {
			layer = "SL"
			layerId = 3
		} else {
			layer = "SO"
			layerId = 2
		}

	} else {
		print "arc3d mismatch at ",secname(),"(",xseg,")"
		layer="ERROR"
		layerId=-1
	}
}

// Do initializations.
// This is repeated for each simulation pass,
// that is, once per synapse location simulated.

proc init() {
	finitialize(Vrest)
	fcurrent()
}


// Run a simulation pass for a single synapse.
// Results of the synaptic stimulation are analyzed
// and summary data is written to a file.

proc runIt() {

	// Set clamps as appropriate for the recording mode.
	// Note that CCinj is not recomputed here since impedance is
	// a complicated procedure and the value should not change
	// from one synapse to the other. Other redundant initializations
	// contained here are to allow command line invocations with
	// different parameters.

	seClamp.dur1 = tError+tInit
	seClamp.amp1 = VCvm

	hold.dur=tError+tInit 
	hold.amp=0	// set to CCinj in runIt()
	hold.del=10	// a short delay for visualization

	if (isCC) {
		seClamp.rs=1e15
		hold.amp=CCinj
	} else {
		seClamp.rs=1
		hold.amp=0
	}

	ns.start=tInit
	tstop=tInit
	cvode.maxstep(10)
	run()

	// Stimulated the synapse and record the
	// results. For this, small time steps
	// are need to get the respone kinetics.
	cvode.maxstep(tDelta) // Set sample rate
	stoprun=0
	continuerun(tInit+tCont) 
	if (stoprun) {
		isError=1
		print "Stopped by request"
		return
	}

 
	// Look at the response and extract peak value,
	// time to peak, and half-height width. If the
	// the initial recording did not capture enough
	// data for the analysis, more run time is allowed
	// in the loop below. First we get time to peak.
	iMin=1e9
	vMax=-1e9
	dMax=-1e9
	tPeak=0
	kPeak=0
	kInit=0
	tRise1=0
	tRise2=0
	tRise=0
	
	// Find the index of the end of the initialization period
	// Note the soma equilibrium values for voltage
	// and current. These are the baselines for EPSP/C.
	for (k=0;k<tvec.size() && tvec.x[k]<=tInit; k+=1) {}
	kInit=k-1
	iRest=ivec.x[kInit]		// rest soma current for VC
	vRest=vvec.x[kInit]		// rest soma voltage for CC
	dRest=dvec.x[kInit]		// rest synapse voltage

	// Make sure that we have passed a peak. If not, keep running
	// until there is an apparent peak or else we have reached the
	// maximum allowed time.
	isDone=0
	while (!isDone) {
		k=tvec.size()-1
		if (k<0) {
			print "Initial tvec is empty at ",secname(),synLoc
			print "This should never occur"
			isError=1
			return
		} else if (isCC==0) {
			for (j=k;!isDone && j>kInit;j=j-1) {
				if (ivec.x[j-1]<=ivec.x[j]) {
					isDone = 1
				}
			}
		} else if (isCC==1) {
			for (j=k;!isDone && j>kInit;j=j-1) {
				if (vvec.x[j-1]<=vvec.x[j]) {
					isDone = 1
				}
			}

		}

		if (isDone) {
			break
		} else if (tvec.x[k]<tInit+tError) {
			print "Continuing simulation at ",secname(),"(",synLoc,")"
			stoprun=0
			continuerun(t+tCont)
			if (stoprun) {
				isError=1
				print "Stopped by request"
				return
			}
		} else {
			print "Failed to detect response peak for ",secname(),synLoc
			isError=1
			return
		}
	}

	// Get peak response voltage at synapse
	for (k=kInit; k<tvec.size(); k+=1) {
		if (dvec.x[k]>dMax) {
			dMax=dvec.x[k]
		}
	}
	
	// For VC, get EPSC peak time, value, and rise time
	if (isCC==0) {
		// Locate the peak (negative relative to rest for EPSC)
		for (k=kInit; k<tvec.size(); k+=1) {
			if (ivec.x[k]<iMin) { // EPSC<0
				iMin=ivec.x[k]
				tPeak=tvec.x[k]
				kPeak=k
			}
		}
		// Get the 20% and 80% or peak rising time values
		for (k=kInit; k<=kPeak; k+=1) {
			if (ivec.x[k]>=0.8*iRest+0.2*iMin) {
				tRise1=tvec.x[k]
				iRise1=ivec.x[k]
			}
			if (ivec.x[k]>=0.2*iRest+0.8*iMin) {
				tRise2=tvec.x[k]
				iRise2=ivec.x[k]
			}
		}
		// Adjust rise time based on actual values at granular time values
		if (iRise2<iRise1) {
			tRise=(tRise2-tRise1)*0.6*(iMin-iRest)/(iRise2-iRise1)
		}
	}

	// For CC, get EPSP peak time, rise time, and value
	if (isCC==1) {
		// Locate the peak (postive relative to rest for EPSP)		 
		for (k=kInit; k<tvec.size(); k+=1) {
			if (vvec.x[k]>vMax) { // Note EPSP>0
				vMax=vvec.x[k]
				tPeak=tvec.x[k]
				kPeak=k
			}
		}
		// Get the 20% and 80% or peak rising time values
		for (k=kInit; k<=kPeak; k+=1) {
			if (vvec.x[k]<=0.8*vRest+0.2*vMax) {
				tRise1=tvec.x[k]
				vRise1=vvec.x[k]
			}
			if (vvec.x[k]<=0.2*vRest+0.8*vMax) {
				tRise2=tvec.x[k]
				vRise2=vvec.x[k]
			}
		}
		// Adjust rise time based on actual values at granular time values
		if (vRise2>vRise1) {
			tRise=(tRise2-tRise1)*0.6*(vMax-vRest)/(vRise2-vRise1)
		}
	}


	if (tPeak==0) {
		isError=1
		print "Failed to detect response peak for ",secname(),synLoc
		print iMin,vMax,tPeak
	}

	// Get the rising (pre-peak) half-height time
	thw1=0
	thw2=0
	isDone=0
	if (isCC==0) for (k=kInit;k<=kPeak && !isDone;k+=1) {
		if (ivec.x[k]<=(iMin+iRest)/2) {
			thw1=tvec.x[k]
			isDone=1
		}
	}
	if (isCC==1) for (k=kInit;k<tvec.size() && !isDone;k+=1) {
		if (vvec.x[k]>=(vMax+vRest)/2) {
			thw1=tvec.x[k]
			isDone=1
		}
	}
	if (!isDone) {
		isError=1
		print "Failed to detect half-rise time for ",secname(),synLoc
	}

	// Now get the falling (post-peak) half-height time.
	// If the current recording does not contain it,
	// extend the simulation for additional time.
	isDone=0
	while (t<tError+tInit && !isDone) {
		if (isCC==0) for (k=kInit;k<tvec.size() && isDone==0;k+=1) {
			if (tvec.x[k]>=tPeak && ivec.x[k]>=(iMin+iRest)/2) {
				thw2=tvec.x[k]
				isDone=1
			}
		}
		if (isCC==1) for (k=kInit;k<tvec.size() && isDone==0;k+=1) {
			if (tvec.x[k]>=tPeak && vvec.x[k]<=(vMax+vRest)/2) {
				thw2=tvec.x[k]
				isDone=1
			}
		}

		// See if we have detected the falling half-height point.
		// If not, run the simulation further to try and locate it.
		if (!isDone) {
			print "Continuing simulation at ",secname(),"(",synLoc,")"
			stoprun=0
			continuerun(t+tCont)
			if (stoprun) {
				isError=1
				print "Stopped by request"
				return
			}
		}
			
	}
	if (!isDone) {
		isError=1
		print "Failed to detect half-fall time for ",secname(),synLoc
	}

	// If error found, do not continue here
	if (isError) return

	// Write the results in summary form to the output file
	// and also print some values on the console log to distract any
	// large semi-aquatic primates watching for their reward cue.
	if (outfile!=nil) {
		outfile.printf("%s,%g",secname(),synLoc)
		outfile.printf(",%s,%g,%g",layer,ory3d,distance(synLoc))
		outfile.printf(",%g,%g",segLen,segArea)
		outfile.printf(",%g",onPPpath)
	}

	if (isCC) {
		print "EPSP (mV) = ",vMax-vRest
		if (outfile!=nil) outfile.printf(",CC,%g",vMax-vRest)
	} else {
		print "EPSC (pA) = ",(iMin-iRest)*1e3
		if (outfile!=nil) outfile.printf(",VC,%g",(iMin-iRest)*1e3)
	}

	print "Rise time (ms) = ",tRise
	print "Peak time (ms) = ",tPeak-tInit
	print "Half-height width (ms) = ",thw2-thw1
	print "Peak synaptic Vm (mV) = ",dMax

	if (outfile!=nil) {
		outfile.printf(",%g",tRise)
		outfile.printf(",%g",tPeak-tInit)
		outfile.printf(",%g",thw2-thw1)
		outfile.printf(",%g,%g",dRest,dMax)
		outfile.printf(",%g,%g,%g",amparNC.weight*1e3,amparSyn.tau1,amparSyn.tau2)
		outfile.printf(",%g,%g,%g",nmdarNC.weight*1e3,nmdarSyn.tau1,nmdarSyn.tau2)
		outfile.printf("\r\n") // Windows (and html) end of line
	}
}


// Run just one synapse in the current section.
proc getSynResp() {local xseg
	xseg = $1
	segLen=L/nseg

	// See where this segment is located
	findLayer(xseg)

	// Skip if this is soma, axon, SL, 
	// or otherwise an error.
	if (layerId<=1 || layerId==3 || layerId>=100) { // SOMA, AXON, SL, Spine etc.
		return
	}

	// Customize the synapse to the layer and apply
	// any randomization rules to synapse properties.
	// Even though Neuron provides normal and log-normal
	// distributions, we do our own conversion from
	// the default distribution for efficiency and
	// in the general interest of knowing exactly what
	// the results will be. Note that exp(x) where
	// x~N(0,s^2) does not have a mean of 1. In general
	// if we want exp(x) to have a mean of 1 we need to
	// have x~N(-0.5*s^2,s^2). In this case the median
	// of exp(x) is not the same as the mean.

	rwm = 1
	rtm = 1
	if (layerId==2 || layerId==4) { // SO and SR
		if (randomizeSyn) {
			rwm = exp(rand.repick()*rcAMPARSynWsd-0.5*rcAMPARSynWsd^2)
			rtm = exp(rand.repick()*rcAMPARSynTsd-0.5*rcAMPARSynTsd^2)
		}
		amparNC.weight=rcAMPARSynWeight*1e-3*rwm
		amparSyn.tau1=minVal(rcAMPARSynTau1,rtm*rcAMPARSynTau2)
		amparSyn.tau2=maxVal(rcAMPARSynTau1,rtm*rcAMPARSynTau2)

		if (randomizeSyn) {
			rwm = exp(rand.repick()*rcNMDARSynWsd-0.5^rcNMDARSynWsd^2)
			rtm = exp(rand.repick()*rcNMDARSynTsd-0.5*rcNMDARSynTsd^2)
		}
		nmdarNC.weight=rcNMDARSynWeight*1e-3*rwm
		nmdarSyn.tau1=minVal(rcNMDARSynTau1,rtm*rcNMDARSynTau2)
		nmdarSyn.tau2=maxVal(rcNMDARSynTau1,rtm*rcNMDARSynTau2)

	} else if (layerId==5) { // LM
		if (randomizeSyn) {
			rwm = exp(rand.repick()*ppAMPARSynWsd-0.5*ppAMPARSynWsd^2)
			rtm = exp(rand.repick()*ppAMPARSynTsd-0.5*ppAMPARSynTsd^2)
		}
		amparNC.weight=ppAMPARSynWeight*1e-3*rwm
		amparSyn.tau1=minVal(ppAMPARSynTau1,rtm*ppAMPARSynTau2)
		amparSyn.tau2=maxVal(ppAMPARSynTau1,rtm*ppAMPARSynTau2)

		if (randomizeSyn) {
			rwm = exp(rand.repick()*ppNMDARSynWsd-0.5*ppNMDARSynWsd^2)
			rtm = exp(rand.repick()*ppNMDARSynTsd-0.5*ppNMDARSynTsd^2)
		}
		nmdarNC.weight=ppNMDARSynWeight*1e-3*rwm
		nmdarSyn.tau1=minVal(ppNMDARSynTau1,rtm*ppNMDARSynTau2)
		nmdarSyn.tau2=maxVal(ppNMDARSynTau1,rtm*ppNMDARSynTau2)

	} else {
		print "Bad layerId value=",layerId
		return
	}

	// Reset synapse conductance and [Mg++] based on
	// types of synapses not blocked. 
	nmdarSyn.mg=lowExtMg			// Assume low [Mg++] as the default
	if (useAMPAR && useNMDAR) {		// Both types of responses together
		nmdarSyn.mg=highExtMg
	} else if (useAMPAR) {			// Only AMPAR
		nmdarNC.weight=0
	} else if (useNMDAR) {			// Only NMDAR
		amparNC.weight=0
	} else {
		print "Both AMPAR and NMDAR are blocked"
		return
	}
 
	// Finish up with the synapses (this is redundant but safety first)
	amparSyn.e=synErev
	nmdarSyn.e=synErev

	// Connect synapses with their postsynaptic element
	// For LM synapses go directly on the dendrite while
	// elsewhere, they are on the spine head if useSpine=1
	synLoc=xseg
	if (useSpine!=1 || layerId==5) {
		amparSyn.loc(synLoc)
		nmdarSyn.loc(synLoc)
	} else {
		// Customize the spine head for the equivalent of
		// a K+ conductance of resistivity gspine.
		// Also, make sure that the spine has an Ra
		// to give a net 500 M-ohm axial resistance.
		spineHead {
			e_pas = (vRest/Rm + gkspine*ekspine)/(1/Rm+gkspine)
			g_pas=1/Rm+gkspine
			Ra=RaSpine
		}
		spineNeck { Ra=RaSpine }

		// Connect synapse to the spine head and then
		// connect the spine neck with the dendrite
		spineHead{ amparSyn.loc(0.5) nmdarSyn.loc(0.5) }
		spineNeck disconnect()
		connect spineNeck(0),synLoc
		soma[somaMidSec] distance(0,0.5)	
	}
	

	// If AOK, go run something useful
	if (!isError) {

		segArea=area(synLoc)
		dvec.record(&v(synLoc))
		iampar.record(&amparSyn.i)
		inmdar.record(&nmdarSyn.i)

		// Indicate where the synapse is located
		printf("\n%s(%g) ",secname(),synLoc)
		printf("layer=%s ory3d=%g dist=%g\n",layer,ory3d,distance(synLoc))
		runIt()

	}

	// If using spines, disconnect the spine and
	// restore distances to the original values.
	if (useSpine==1) {
		spineNeck disconnect()
		soma[somaMidSec] distance(0,0.5)
	}
}

// Simulate synaptic responses with a synapse
// placed in each segment, one at a time for all
// segments in SO, SR, or LM layers.

proc sampleSynResps() {

  outfile = new File()
  outfile.wopen(outPath)

  // Write a header line to outfile
  outfile.printf("%s","section,loc,layer,y,dist,len,area,trunk,type")
  outfile.printf("%s",",peakValue,riseTime,peakTime,halfWidth,synRest,synPeak")
  outfile.printf("%s",",gAMPAR,tau1AMPAR,tau2AMPAR,gNMDAR,tau1NMDAR,tau2NMDAR")
  outfile.printf("\r\n") // Windows (and html) end of line

  // Loop through all sections
  forall {

	// Loop through segments placing the synapse in
	// appropriate segments and recording the synaptic
	// response at that location.

	for (xseg,0) {

		// Get (and write) the response for a synapse
		// at location xseg in the current section.
		getSynResp(xseg)
	}
  }

  outfile.close()
}

// Simulate a single synaptic activation and write the results
// to a file. This function would normally only be invoked manually.
// Invocation is saveSynResp(xseg). The path name of the file to save
// to is taken from global string savePath. Voltages are saved as mV
// and currents as pA. Time is in units of ms.

proc saveSynResp() {local xseg
	xseg=$1
	saveFile = new File()
	saveFile.wopen(savePath)

	tCont=200 // make sure whole response is captured
	getSynResp(xseg)
	
	// Write a CSV header line
	saveFile.printf("time,soma,dend,iampar,inmdar\r\n")

	// Write the appropriate data values.
	if (isCC) { // save time and voltage, etc.
		for(k=0; k<tvec.size(); k+=1) {
			saveFile.printf("%g,%g,%g",tvec.x[k],vvec.x[k],dvec.x[k])
			saveFile.printf(",%g,%g\r\n",1e3*iampar.x[k],1e3*inmdar.x[k])
		}
	} else { // save time and current, etc.
		for(k=0; k<tvec.size(); k+=1) {
			saveFile.printf("%g,%g,%g",tvec.x[k],1e3*ivec.x[k],dvec.x[k])
			saveFile.printf(",%g,%g\r\n",1e3*iampar.x[k],1e3*inmdar.x[k])
		}
	}

	saveFile.close()
}


// ================================
// End Subroutines
// ================================


// Find out how many sections make up the soma.
// Then look for the one with the largest diameter
// and designate it as the middle section.
forsec "soma" {nSomaSec=nSomaSec+1}
maxDiam=0
for (k=0;k<nSomaSec;k=k+1) soma[k] {
	if (maxDiam<=diam) {
		maxDiam=diam
		somaMidSec=k
	}
}
print "nSomaSec=",nSomaSec," somaMidSec=",somaMidSec	 

// In Neuromorpho SWC files, the soma is sometimes
// represented by a single point with the radius being
// the radius of a sphere of equivalent area. When
// CVAPP converts these files to HOC, the single point
// becomes two points with a small (0.01) difference in
// Z coordinates. Detect this problem and change the
// soma to a cylinder with area equivalent to the
// sphere intended in the SWC file. The cylinder is
// arbitrarily aligned along the Y axis.

access soma[somaMidSec]

if (nSomaSec==1) {
if (n3d()==2) {
if (abs(z3d(0)-z3d(1))<0.02) {
if (diam3d(0)==diam3d(1)) {
if (y3d(0)==y3d(1)) {
if (x3d(0)==x3d(1)) {
	r = diam3d(0)/2
	pt3dchange(0,x3d(0),y3d(0)-r,z3d(0),2*r)
	pt3dchange(1,x3d(1),y3d(1)+r,z3d(1),2*r)
	print "Soma area changed to ",4*PI*r^2
}}}}}}


// Set nseg so that each segment is at most seglen long.
// For this purpose, an odd number of segments is not
// important since we are not going to change nseg again.
// soma area(0.5) // from fixnseg.hoc -- what does it do?
forall {
	nseg = int(1-1e-8+L/maxSegLen)
}



// Initialize the lookup table for detecting apical_dendrites
// that connect somewhere into PP. Each section in LM is found
// and that section and all parents in turn are flagged in adOnPPpath.
// This assumes that the most distal part of a section is in LM
// if any part of the section is (this should always be true).
forsec "apical_dendrite" {
	findLayer(1) // test most distal part of the section
	if (layerId==5) {
		adOnPPpath.x(curSecnum)=1			
		while (curSecnum>=0) {
			adOnPPpath.x(curSecnum)=1
			apical_dendrite[curSecnum] {
				secRef = new SectionRef()
				secRef.parent parseSecname()
			}
			if (strcmp(curSecname,"apical_dendrite")!=0) break
		}
	}	
}

// Create any special mechanisms. The soma is the
// context for the clamps, but the synapse is
// later relocated to different places as needed.

soma[somaMidSec] {

	// Set up a voltage clamp to mimic experimental
	// conditions prior to stimulation.
	seClamp = new SEClamp(0.5)
	seClamp.dur1 = tError+tInit
	seClamp.amp1 = VCvm

	// Just for fun, benchmark how long it takes for VC
	// to propagate from the soma into dendrites.
	// seClamp.dur1 = 500
	// seClamp.amp1 = -65
	// seClamp.dur2 = 500
	// seClamp.amp2 = -80

	hold = new IClamp(0.5)
	hold.dur=tError+tInit 
	hold.amp=0	// set to CCinj in runIt()
	hold.del=10	// a short delay for visualization

	// Adjust for (initial) recording mode
	// This starts things out right when
	// running in interactive mode. 	
	if (isCC) {
		seClamp.rs=1e15
	} else {
		seClamp.rs=1
		hold.amp=0
	}

	// Initialize the stimulation and synapse
	ns = new NetStim(.5)
	ns.start=0 			// start time will be reset later
	ns.number=1			// only a single spike is simulated
	ns.interval=1e9 		// not really used here

	// Create a synapse that will be moved around.
	// Even though AMPAR and NMDAR actually coexist
	// in one synapse, for the present purpose they
	// are treated separately. Note that this is not
	// entirely valid if a spine neck is involved
	// since the membrane potential at the spine
	// head, where the receptors are located, is not
	// the same as on the main dendrite membrane. Hence
	// AMPAR induced depolarizations are off a bit.
	// However, the parameters to accurately model
	// the spine are not well determined experimentally. 
	amparSyn = new Exp2Syn(0.5)
	nmdarSyn = new Exp2NMDAR(0.5)	
	amparSyn.e=synErev
	nmdarSyn.e=synErev
}
soma[somaMidSec] distance(0,0.5)	// sets soma midpoint as origin


// Set up to record from the soma
ivec.record(&seClamp.i)
vvec.record(&soma[somaMidSec].v(0.5))
tvec.record(&t)

// Create network connections. These are customized
// based on location in getSynResp().
amparNC = new NetCon(ns, amparSyn, 0, 0, 0)
nmdarNC = new NetCon(ns, nmdarSyn, 0, 0, 0)


// Insert normal mechanisms including any active channels
if (addChannels && useHemond2008) { // Mimic the Hemond et al. model

	// Add channels to different types of sections
	// For conducances known to be 0, channels are not added.
	// Note that gKc and gahp are both dependent on internal
	// calcium concentration. The use of cacum here is from
	// the Hemond et al. 2008 model and is used AS IS. It is
	// not intended to be highly accurate with respect to
	// the underlying biophysics of calcium dynamics.
	
	forsec "dendrite" { // both apical and basal dendrites
		if (gna>0) {insert na3 gbar_na3=gna sh_na3=sh}
		if (ghd>0) {insert hd  ghdbar_hd=ghd sh_hd=sh}
		if (gc>0 || gKc>0) {insert cacum depth_cacum=diam/2}
		if (gc>0) {
			insert cal gcalbar_cal=gc
			insert can gcanbar_can=gc
			insert cat gcatbar_cat=gc
		}
		if (gkdr>0) {insert kdr gkdrbar_kdr=gkdr sh_kdr=sh}
		if (gkap>0) {insert kap gkabar_kap=gkap sh_kap=sh}
		if (gKc>0)  {insert cagk gbar_cagk=gKc}
		if (gahp>0) {insert KahpM95 gbar_KahpM95=gahp}
	}
	forsec "soma" {
		insert na3 gbar_na3=gna sh_na3=sh
		insert hd  ghdbar_hd=ghd sh_hd=sh
		insert cacum depth_cacum=diam/2
		insert cal gcalbar_cal=gc
		insert can gcanbar_can=gc
		insert cat gcatbar_cat=gc
		insert kdr gkdrbar_kdr=gkdr sh_kdr=sh
		insert kap gkabar_kap=gkap sh_kap=sh
		insert km gbar_km=gkm sh_km=sh
		insert kd gkdbar_kd=gkd sh_kd=sh
		insert cagk gbar_cagk=gKc
		insert KahpM95 gbar_KahpM95=gahp

	}
	forsec "axon" {
		insert na3 gbar_na3=gnaAxon sh_na3=sh
		insert kdr gkdrbar_kdr=gkdr sh_kdr=sh
		insert kap gkabar_kap=gkap sh_kap=0 // Note 0 mV shift
	}
}

if (addChannels && !useHemond2008) { // Only Kdr and K-A with no shifts

	forsec "dendrite" { // both apical and basal dendrites
		insert kdr gkdrbar_kdr=0 sh_kdr=0
		insert kap gkabar_kap=0
	}
	forsec "soma" {
		insert kdr gkdrbar_kdr=gkdr sh_kdr=0
		insert kap gkabar_kap=gkap sh_kap=0

	}
	forsec "axon" {
		insert kdr gkdrbar_kdr=gkdr sh_kdr=0
		insert kap gkabar_kap=gkap sh_kap=0
	}
}


// Customize each compartment as needed
forall {

	insert pas	// Add passive conductances

	// Set passive params section variables. Ra is special and
	// unlike pas properties, cannot change within a section.
	Ra=RaAll

	// Set reversal potentials based on what channels are found
	if (ismembrane("na3")) {
		ena=55		// Na+ reversal
	}
	if (ismembrane("kdr") || ismembrane("kap") || ismembrane("kad")) {
		ek=-90		// K+ reversal
	}

	// Handle segments one at a time to allow for any
	// distance or layer dependencies.
	for (xx,0) {

		findLayer(xx)

		// Accumulate membrane area based on layerId
		if (layerId==0) { 
			totSoma += area(xx)
		} else if (layerId==2) { 
			totBasal += area(xx)
		} else if (layerId>=3) {
			totApical += area(xx)
		} 

		// Get spine adjustment factors for this location
		spineAdj=1
		if (layerId==5) { spineAdj=PPSpineAdj }
		if (layerId==2 ||layerId==4) { spineAdj=RCSpineAdj }

		xxDist=abs(distance(xx))

		// Start with baseline values
		cmHere = Cm
		gmHere = 1/Rm

		if (layerId>1) { // not soma or axon

			// Adjust membrane capacitance
			cmHere *= spineAdj

			// Adjust membrane conductance
			gmHere *= spineAdj
		}

		if (layerId==1) { // Axon
			raHere = RaAxon
		}

		// Set location dependent passive params
		e_pas(xx) = Vrest
		cm(xx) = cmHere
		g_pas(xx) = gmHere		

		// Set any active channel conductances that are dependent
		// on layer or distance from soma. This is not used for
		// the Hemond et al. 2008 model.
		if (addChannels && !useHemond2008) {
			if (layerId==2 || layerId==4 || layerId==5) { // SO or SR or LM
				if (ismembrane("kap")) {
					gkabar_kap(xx) = gkap*(1+gkaSlope*minVal(xx,gkaSlopeMaxDist))
				}
				if (ismembrane("kad")) {
					gkabar_kad(xx) = gkap*(1+gkaSlope*minVal(xx,gkaSlopeMaxDist))
				}
				gkdrbar_kdr(xx) = gkdr
			}
		}
	}
}


// If CC, get the input impedance of the cell
// so that we can inject the right amount of current
// to get to the correct holding potential.


if (isCC) {
	init() // Set initial voltages
	impedance=new Impedance(0)

	soma[somaMidSec] {
		impedance.loc(0.5)
//		impedance.compute(0,1) // this fails in Neuron 7.1
		impedance.compute(0) // good enough for now
	}

	rinput=impedance.input(0)
	print "Input impedance = ",rinput
	CCinj=(CCvm-soma[somaMidSec].v(0.5))/rinput
}

// Gather some morphology statistics
forall for (xseg,0) {
	findLayer(xseg)
	segLen=L/nseg
	if (layerId==0) {
		totSomaLen += segLen
		totSomaArea += area(xseg)
	}
	if (layerId==1) {
		totAxonLen += segLen
		totAxonArea += area(xseg)
	}
	if (layerId==2) {
		totSOlen += segLen
		totSOarea += area(xseg)
	}
	if (layerId==4) {
		totSRlen += segLen
		totSRarea += area(xseg)
		if (onPPpath) {
			totSRTlen += segLen
			totSRTarea += area(xseg)
		}
	}
	if (layerId==5) {
		totLMlen += segLen
		totLMarea += area(xseg)
	}
}

// Initialize the random number generator with a seed.
// Even though this is not needed for Random, explicitly
// start with a normal distribution N(0,1).
// The total membrane area is included in the seed
// to allow each cell some distinctiveness.
rand = new Random(randSeed+totSomaArea+totLMarea+totSRarea+totSOarea)
rand.normal(0,1)


// Run the simulation if requested.
if (runStim) {
	sampleSynResps()
}

// Remind us again what we were running when all this
// started (in some cases, a long long time ago).
print " "
print "Cell = ",cell
print "Out = ",outPath
print "Maximum segment length = ",maxSegLen
if (useAMPAR) print "AMPAR responses are simulated"
if (useNMDAR) print "NMDAR responses are simulated"

// Show morphology statistics
print " "
print "Total Soma length and area = ",totSomaLen," ",totSomaArea
print "Total Axon length and area = ",totAxonLen," ",totAxonArea
print "Total LM length and area = ",totLMlen," ",totLMarea
print "Total SR length and area = ",totSRlen," ",totSRarea
print "Total SR trunk length and area = ",totSRTlen," ",totSRTarea
print "Total SR non-trunk length and area = ",totSRlen-totSRTlen," ",totSRarea-totSRTarea
print "Total SO length and area = ",totSOlen," ",totSOarea
print " "