/* Exploring the effects of dendritic morphology on spike time patterns. */
/* The goal of this exercise is both to explore Mainen and Sejnowski's model
(for more exploration ideas see the hands-on exercise modeldb.html)
and to investigate the transition from quiescence to bursting using injected
DC current as a bifurcation parameter. */

/* Instructions:

1. Read the abstract of Mainen & Sejnowski's 1996 Nature paper.  (The
PDF is on the course website.)

2. Download and compile the 1996 Mainen & Sejnowski model from ModelDB.

3. Run mosinit.hoc and explore the behavior of the different cell types.

4. Quit NEURON, restart mosinit.hoc afresh, and select a cell type (e.g.
L3 Aspiny)

Then run the following code (e.g. cut and paste the line below into the
hoc window)

xopen("...pathname to wherever you put the file mainen-sejnowski-1996-problem.hoc...")

For me, this would translate as:

xopen("/Users/pjthomas/Documents/Courses/CNS/CaseCNS/b378s08/b378s08hwk/hwk3/mainen-sejnowski-1996-problem.hoc")

5. For the other cell types, modify the hoc code to change the range of
currents used to span from zero to 1.5 times the default current used in
the original model.

*/

cvode_active(0) // This shuts off variable timestep.  Alternatively, by hand you can 
// check Main Menu / Tools / Variable Time Step to make sure variable dt is NOT checked.

// make a list of different current injection amplitudes to use
objref IClamp_amps,freq, gplot
IClamp_amps = new Vector()

/* These are examples of short cuts that might be useful for debugging later...
//IClamp_amps.append(0,) // speed things up for debugging -- just do a single amplitude = 0 
//IClamp_amps.append(0,.005,.010) // speed things up for debugging -- three amplitudes, no spikes 
//IClamp_amps.append(0,.025,.050) // speed things up for debugging -- three amplitudes, some with spikes 
//tstop=100			  // speed things up for debugging -- go to t=100 instead of 1000
*/

//// For L3 Aspiny (default drive is 0.05) -- v
//IClamp_amps.indgen(0,.075,.005) // vec.indgen(vmin,vmax,dv) -- same as "0:.005:.075" in matlab
IClamp_amps.append(0,.025,.050) // speed things up for debugging -- three amplitudes, some with spikes 
tstop=200

// for each value of IClamp_amps, run a simulation and keep track of the asymptotic frequency
freq = new Vector()
freq.resize(IClamp_amps.size)

objref vec, bin, times, isi, freqtmp, ones
vec = new Vector()		// Vector to hold the voltage
bin = new Vector()		// 1 if a spike occurred in a given time step, else 0
times = new Vector()	// spike times
isi = new Vector()		// Inter spike intervals (ISI)
freqtmp = new Vector()		// local spike "frequency" i.e. 1/ISI
ones = new Vector()		// a vector of all ones, the size of isi.

for ii_amp=0,IClamp_amps.size-1 {
	print "Main Loop Top: ii_amp = ",ii_amp
	IClamp[0].amp=IClamp_amps.x[ii_amp] // set IClamp amplitude for this run
	vec.record(&soma.v(.5),dt) // prepares vector "vec" to record soma.v(.5) during the simulation
	init()					//
	run()					// 
	bin.spikebin(vec,0) 	// dest_vec.spikebin(data,threshold (mv)) finds the threshold crossings
	times.indvwhere(bin,"==",1) // find those indices at which a spike occurred.
	//print "starting: 	if (times.size < 2) {	, etc."
	// if < 2 spikes then call frequency 0, otherwise calculate ISIs and frequency
	if (times.size < 2) {print "no ISIs"
		freq.x[ii_amp] = 0 } else {		
		print "at least one ISI"
		times.mul(dt) 			// multiply each element of times by dt to give the actual spike times
		isi.resize(times.size-1)// expand ISI vector to appropriate size 
		for ii=0,times.size-2 {	isi.x[ii]=(times.x[ii+1]-times.x[ii]) }
		ones.resize(isi.size)	// create a vector the same size as isi
		ones.fill(1,0,isi.size-1) // create a vector the same size as isi and fill with "1"s.
		freqtmp = ones.div(isi)	// find local frequency as 1/isi
		freq.x[ii_amp] = freqtmp.x[freqtmp.size-1]} // take last frequency value as "asymptotic" -- crude but it's a start
	print "Main Loop Bottom: ii_amp = ",ii_amp," freq = ",freq.x[ii_amp] } // end of "for ii_amp=0,IClamp_amps.size"
print "main loop completed"

// Print the results to the Terminal window.
for ii=0,freq.size-1 {print IClamp_amps.x[ii],freq.x[ii]}
// Graph the results
objref g
g = new Graph()
g.size(IClamp_amps.x[0],IClamp_amps.x[IClamp_amps.size-1],-.1*freq.max,1.1*freq.max)
freq.plot(g,IClamp_amps,3,5)

/* Debugging note (applies to mac, maybe also to linux/unix):

On the mac, you drag and drop a hoc file (say, mosinit.hoc) onto the
nrngui.app icon to start NEURON. When NEURON is done -- including if it
crashes -- the Terminal window it was running in disappears by default. 
This prevents seeing the final error message that preceeded the crash.

To get around this, you can copy the following command into the Terminal
window:

bash$ cd ; clear; cp .nrn_as2sh /tmp/nrn_as2sh$$; sh /tmp/nrn_as2sh$$;rm -f /tmp/nrn_as2sh$$

(don't type "bash$", that's just the prompt).  The file .nrn_as2sh in
your home directory might contain something like this:

#!/bin/sh
NRNHOME="/Applications/NEURON-6.1/nrn"
NEURONHOME="${NRNHOME}/share/nrn"
NRNBIN="${NRNHOME}/umac/bin/"
PATH="${NRNHOME}/umac/bin:${PATH}"
export NRNHOME
export NEURONHOME
export NRNBIN
export PATH
nrncarbon=yes
export nrncarbon
cd "${NRNHOME}/umac/bin"
./nrngui.sh "/Users/you/whereveryouputthemodfiles/patdemo/mosinit.hoc"

This file is recreated each time you drag and drop mosinit.hoc (or
another file) onto the NEURON gui icon.  If you call NEURON directly
from the Terminal as above, the Terminal window stays around after you
crash NEURON while debugging, which makes it easier to find the errors.
 
 */