/*************************************
*
*	Modified from fitFR.hoc in the Coskren/July2012/baseline directory,
*	to apply Aniruddha's conductances and axon to all of Patrick's neurons.
*
*************************************/

objref voltage_vec, time_vec, dendritic
load_file("fixnseg.hoc")
load_file("readcell_nomechanisms.hoc")

// *** Globals ***

// Reversal potential for the 'pas' membrane mechanism.  NEURON defaults to
// -70 mV for this; since we want to assume standard membrane properties for
// all neurons (in order to isolate differences due to morphology), this is
// fine for our purposes.
E_PAS = -80	// Aniruddha's value
PHI_DFLT = 52 / 2e-3
BETA_SOMA = 1/100
BETA_DEND = 1/20
CEILING_CA = 1000
AX_GNASCALE = 3		// How do axonal gNa values compare to soma?

// Neuron uses nA for IClamp; we generally use pA
kPicoToNanoMultiplier = 0.001

// The soma created here will be replaced by read_cell, but it's needed in
// order for some of the following code to be accepted by the Hoc interpreter.
create soma

flag_spines = 1

objref voltageClamp
objref stim, s

// Global variables from Christina's aux_procs.hoc script.
objref FRout
strdef refStr
objref exptVec, iVec
objref ihold
objref tree_root
// copied from coskren_procs.hoc
STD_SOMA_RADIUS = 8.28  // Only used for electrotonic measurements
STD_SOMA = STD_SOMA_RADIUS
CM = 0.83382966	// used by Aniruddha
RA = 150
celsius = 37	// used by Aniruddha; Coskren used 21

// Global variables from Christina's init_model.hoc script.
spinescale = 1
dendscale = 1

// Global variables from Christina's rec_volt_justV.hoc script
strdef fname
objref tVec, vs, caVec, iVec
objref vs0, vd0, vd1, vd2, vd3	// voltages at the connecting ends of each section
objref sref
objref camVec, cacVec, catVec		// [ca] in middle shell, core, & total
objref caB					// buffered Ca in outer shell
objref aVec, xVec			// for ATPase & exch vec
objref kdrV, nafV, kahpV, cahvaV
objref kslowV, kaV, hV, napV, calvaV
objref canV				// for CAN current
objref spiketimes
objref apc, isi, fr				// APcount
objref fout
objref stVec			// vector containing applied current steps
objref alCAN			// alpha channel opening for CAN
objref icangraph, ecaV, ipmpV

// Global variables from Christina's custominit.hoc script

// This is the current in nA that must be applied at the injection site
// in order to hold v at that location to the desired potential.
IHOLD = 0

// Global variables from Christina's batchrun.hoc script
objref ihold // an IClamp used to deliver the holding current

// Global variables from Christina's mainMac_PFC_wSEClamp.hoc script
INITDUR = 80
V0 = -70

// Global variables containing neuron locations.  The variables neuron_paths
// and neuron_names are parallel lists of strings, such that for any index i,
// neuron_paths[i] is the full path to the neuron, and neuron_names[i] is its
// name.  (Technically, the values are, for example, neuron_paths.o(i).s(),
// since Hoc List objects have no [] operator and these lists contain String
// objects rather than raw Hoc strings.)
objref neuron_paths
objref neuron_names
objref neuron_resting_potentials
neuron_paths = new List()
neuron_names = new List()
objref neuron_path_str
objref neuron_name_str

/*******************************************************
 * Establish neuron locations.  (This is ugly, but I don't want to try to set
 * up reading them from a file in Hoc.)
 */
/* Begin old neurons */

neuron_path_str = new String("models/Aug3CellA-all-apicalbasal-spiny.hoc")
neuron_name_str = new String("Aug3CellA")
neuron_paths.append(neuron_path_str)
neuron_names.append(neuron_name_str)

neuron_path_str = new String("models/Aug3_Slice1CellC-all-apicalbasal-spiny.hoc")
neuron_name_str = new String("Aug3_Slice1CellC")
neuron_paths.append(neuron_path_str)
neuron_names.append(neuron_name_str)

neuron_path_str = new String("models/Aug3_2006CellE-all-apicalbasal-spiny.hoc")
neuron_name_str = new String("Aug3_2006CellE")
neuron_paths.append(neuron_path_str)
neuron_names.append(neuron_name_str)

neuron_path_str = new String("models/Aug3_2006CellF-all-apicalbasal-spiny.hoc")
neuron_name_str = new String("Aug3_2006CellF")
neuron_paths.append(neuron_path_str)
neuron_names.append(neuron_name_str)

neuron_path_str = new String("models/Aug3_2006CellG-all-apicalbasal-spiny.hoc")
neuron_name_str = new String("Aug3_2006CellG")
neuron_paths.append(neuron_path_str)
neuron_names.append(neuron_name_str)

neuron_path_str = new String("models/Feb27_IR2n-all-apicalbasal-spiny.hoc")
neuron_name_str = new String("Feb27_IR2n")
neuron_paths.append(neuron_path_str)
neuron_names.append(neuron_name_str)

/* End old neurons */

/* Begin young neurons */

neuron_path_str = new String("models/Dec15_2006CellE-all-apicalbasal-spiny.hoc")
neuron_name_str = new String("Dec15_2006CellE")
neuron_paths.append(neuron_path_str)
neuron_names.append(neuron_name_str)

neuron_path_str = new String("models/Jun7_IR1d-all-apicalbasal-spiny.hoc")
neuron_name_str = new String("Jun7_IR1d")
neuron_paths.append(neuron_path_str)
neuron_names.append(neuron_name_str)

neuron_path_str = new String("models/May3_IR2d-all-apicalbasal-spiny.hoc")
neuron_name_str = new String("May3_IR2d")
neuron_paths.append(neuron_path_str)
neuron_names.append(neuron_name_str)

neuron_path_str = new String("models/May3_IR2h-all-apicalbasal-spiny.hoc")
neuron_name_str = new String("May3_IR2h")
neuron_paths.append(neuron_path_str)
neuron_names.append(neuron_name_str)

neuron_path_str = new String("models/May3_IR2i-all-apicalbasal-spiny.hoc")
neuron_name_str = new String("May3_IR2i")
neuron_paths.append(neuron_path_str)
neuron_names.append(neuron_name_str)

neuron_path_str = new String("models/May3_IR2t-all-apicalbasal-spiny.hoc")
neuron_name_str = new String("May3_IR2t")
neuron_paths.append(neuron_path_str)
neuron_names.append(neuron_name_str)

/* End young neurons */
/* Done establishing neuron locations. ********************/

// *** Logging ***
// Just a handy function and some globals to make it possible to turn logging on
// and off.
LOG_NONE = 0
LOG_ERROR = 1
LOG_WARNING = 2
LOG_INFO = 3
LOG_TRACE = 4

LOG_LEVEL = 1

// Logs the string passed to it iff the global LOG_LEVEL is >= the log level
// provided as an argument.
//   The log string is passed to printf, so should contain an explicit \n
// character if a newline is desired.
// $1: the log level to compare to LOG_LEVEL
// $s2: the string to print
proc Log() { local log_level localobj log_string
  log_level = $1
  log_string = new String($s2)

  if (LOG_LEVEL >= log_level) {
    printf(log_string.s())
  }
}

/**************************************************
 * The following functions are based on the contents of Christina's
 * add_axon.hoc script.
 */

proc add_yadav_axon() {
  Log(LOG_TRACE, "Entering method add_yadav_axon\n")

  xopen("aniruddha_young10axon.hoc")
  Axon = 1
  define_shape()

  // note that channels have not been added.

  Log(LOG_TRACE, "Leaving method add_yadav_axon\n")
}


/**************************************************
 * The following functions are based on Aniruddha's 
 * model_setup*.hoc and young_start_act_5params.hoc files.
 */


proc define_pass() {

    set_passive($1)
    init()
}

/*****************************************************

    set_passive

    input	$1	1 or 0:  use pasR, not pas?

*****************************************************/
proc    set_passive() {

    geom_nseg(100,0.1)
    v_init  = V0

    forall {
	// set passive params
	Ra = RA
	cm = CM	

        insert pas 
        // reversal potentials
	e_pas = E_PAS
    }
}
	

proc init_yadav_model() {
  Log(LOG_TRACE, "Entering method init_yadav_model\n")

        define_pass(0)
	set_passive(0)
	make_active_dendrite()
	linear_migliore_ratios(1, -0.2, -0.2)

	adj_ka(0.002)
	adj_k2(0.0001)
	adj_h(HVAL)	// not used originally by Aniruddha, but fitted in later simulations.  Default 0.0001
	adj_caL(0.00025)
	adj_nap(0.000003)
	adj_caT(0)

	adj_kca(0.25)
	adj_km(0.0051)
	adj_kahp(0.0001)
	adj_skahp(0.004)
	adj_caL(5.41105e-05)
	adj_naf(0.11032)
	adj_kdr(0.08)

	forall {ek=-95}
	forall {ena=50}
	forall {vrev_naf=-3.5}
	forall {vrev_kdr=29.5}

	forall {if ( ismembrane("cad")) phi_cad=5000}

  Log(LOG_TRACE, "Leaving method init_yadav_model\n")
}



/**************************************************
 * The following functions are based on the contents of Christina's
 * rec_volt_justV.hoc script.
 */

/*
 * Records a vector of AP times (apc) at the center of the currently accessed
 * section, as well as the discrete simulation times (tVec) and soma voltage
 * (vs).  The variables apc, tVec and vs are all global.
 */
proc set_dataVec() {
  Log(LOG_TRACE, "Entering method set_dataVec\n")
  spiketimes = new Vector()
	apc = new APCount(0.5)
	apc.record(spiketimes)
  vs = new Vector()  // Voltage at soma
	tVec = new Vector()  // Vector of time steps
	tVec.record(&t)
	vs.record(&soma.v(0.5))
  Log(LOG_TRACE, "Leaving method set_dataVec\n")
}

/*************************************************************
  volt2txt

  Write time and somatic voltage vectors to a human-readable file named
  '$s1_Vonly.txt'.

  Arguments:
    $s1	basename of output file '$s1_Vonly.txt'

  The following data are written to the output file:

  If stim exists, its amp, del, dur
  time
  soma voltage
*************************************************************/
proc volt2txt() {
  Log(LOG_TRACE, "Entering method volt2txt\n")
	zero = 0
	fout = new File()
	sprint(fname, "%s_Vonly.txt", $s1)
  fout.wopen(fname)
//  printf("volt2bin in rec_volt_simple.hoc:  Opened fname ->%s<-, with $s1 = ->%s<-\n", \
//         fname, $s1)
  if(name_declared("stim")) {
    x = 1  //number of pulses
    fout.printf("Voltage trace.  #pulses: %d amp:%f del:%f dur:%f\n", x, \
                stim.amp, stim.del, stim.dur)
    }
	// set up all vectors for reading
	set_dataVec()

	init()
	Log(LOG_TRACE, "Leaving method run()\n")
	run()
	Log(LOG_TRACE, "Leaving method run()\n")

  // "vs" -> voltage at soma.  Global, defined in set_dataVec.
  fout.printf("%d  # Vector size\n", vs.size())
	vs.printf(fout)

  fout.close()
  Log(LOG_TRACE, "Leaving method volt2txt\n")
}

/*************************************************************
  volt_nobin

  Set up all the interesting vectors as in volt2txt(), but do
  not write the binary file.
*************************************************************/

proc volt_nobin() {
  Log(LOG_TRACE, "Entering method volt_nobin\n")
	zero = 0

	// set up all vectors for reading
	set_dataVec()

	init()
	run()
	if(apc.n > 1) {
	  calcFR_bounds(0,tstop)
  }
  Log(LOG_TRACE, "Leaving method volt_nobin\n")
}

/********************************************************
  calcFR_bounds()

  calculate the mean FR and CV during a specified time
  window.

  Globals:
    apc: "Action Potential Count", created in setDataVec.

  Arguments:
    float $1  left endpoint of time window
    float $2  right endpoint of time window
********************************************************/
proc calcFR_bounds() { local k, tmx
  Log(LOG_TRACE, "Entering method calcFR_bounds\n")
  objref isi, fr  // Globals

  isi = new Vector()
  fr = new Vector()

  // printf("APC: %f\n", apc.n)
  // printf("Spiketimes size: %f\n", spiketimes.size())

  for(k = 0; k < apc.n - 1; k = k + 1) {
    if( spiketimes.x[k] >= $1 && spiketimes.x[k+1] <= $2) {
      isi.append(spiketimes.x[k + 1] - spiketimes.x[k])
      fr.append(1000 / isi.x[isi.size - 1])
    }
  }
  if(fr.size == 0) {
//    printf("Found %d spikes; FR mean = 0, stdev 0, CV 0\n", apc.n)
    return
  }
  if(fr.size > 2) {
//    print "FR mean = ", fr.mean, " stdev ", fr.stdev, " CV ", fr.stdev / fr.mean
  } else {
//    printf("Found %d spikes; FR mean = %.1f\n", apc.n, fr.mean)
  }
  Log(LOG_TRACE, "Leaving method calcFR_bounds\n")
}

/******************************************************************************
  eval_FRandCV()

  Evaluates the firing rate and coefficient of variation for a time series.

  Globals:
    stim must be defined

  Inputs:`
	  float  $1	start time for FR / CV window
    float  $2	end time for FR / CV window
    float  $3	amplitude of current injection
    strdef $s4 file basename
    int    $5	0 or 1, write .Vbin file?

*******************************************************************************/
func eval_FRandCV() {  local old_tstop, old_dur, mnFR, old_del
  Log(LOG_TRACE, "Entering method eval_FRandCV\n")
  old_tstop = tstop
  old_dur = stim.dur

  tstop = $2
  stim.amp = $3
  if(stim.dur == 0) {
    stim.dur = tstop
  }

  if($5) {
    volt2txt($s4)
    Log(LOG_TRACE, "\tDone volt2bin()\n")
  } else {
    volt_nobin()
  }

  calcFR_bounds($1, $2)
  if($5) {
    Log(LOG_TRACE, "\tDone calcFR_bounds() \n")
  }

  tstop  = old_tstop
  stim.dur = old_dur

  if(fr.size > 0) {
    mnFR = fr.mean
  } else {
    mnFR = 0
  }

  return mnFR
  Log(LOG_TRACE, "Leaving method eval_FRandCV\n")
}

/******** End functions from rec_volt_justV.hoc ***********/

/**************************************************
 * The following functions are based on the contents of Christina's
 * aux_procs.hoc script
 */

proc set_epas_negative() {
  Log(LOG_TRACE, "Entering method set_epas_negative\n")
  forall e_pas = -1 * $1
  Log(LOG_TRACE, "Leaving method set_epas_negative\n")
}

// $1: Half-activation voltage shift for Na kinetics
proc shift_NaKin() {
  Log(LOG_TRACE, "Entering method shift_NaKin\n")
  forall if(ismembrane("na_ion"))	vshift_na = $1
  Log(LOG_TRACE, "Leaving method shift_NaKin\n")
}

// $1: Half-activation voltage shift for Na kinetics
proc shift_NaSlopes() {
  Log(LOG_TRACE, "Entering method shift_NaSlopes\n")
  forall if(ismembrane("na_ion"))	{
	qa_na = $1 * 9
	qi_na = $1 * 5
  }
  forall if(ismembrane("kv"))	{
	qa_kv = $1 * 9
  }
  Log(LOG_TRACE, "Leaving method shift_NaSlopes\n")
}
/******** End functions from aux_procs.hoc ***********/


objref axonal, dendritic, somadendrite, apical, basal

proc setup_SecLists() {

    axonal = new SectionList()
    dendritic = new SectionList()
    somadendrite = new SectionList()
    apical = new SectionList()
    basal = new SectionList()

    forsec "soma" {
        somadendrite.append()
    }
    forsec "apical" {
        dendritic.append()
        apical.append()
        somadendrite.append()
    }
    forsec "basal" {
        dendritic.append()
        basal.append()
        somadendrite.append()
    }
    forsec "axon" {
        axonal.append()
    }
    forsec "AxonInitseg" {
        //axonal.append()
        dendritic.append()  	// Axon Init Seg was added to 'dendritic' in Aniruddha's setup_model*.hoc file.
    }
}


/**************************************************
 * The following function is based on the contents of Christina's
 * custominit.hoc script
 */

proc init() { local dtsav, tstopsav, temp
  Log(LOG_TRACE, "Entering method init()\n")
  finitialize(v_init)
  dtsav = dt
  dt = 0.05  // or something larger if stability and accuracy are OK
  t = -1e4
  tstopsav = tstop
  tstop = t + INITDUR
  temp = cvode.active()
  if (temp != 0) {
    cvode.active(0)
  }

  /*****
  voltageClamp.rs = 0.01
  voltageClamp.toff = 0
  voltageClamp.amp = V0
  *****/

  while (t < tstop) {
    fadvance()
  }

  /****
  IHOLD = voltageClamp.i
  printf("In custom init.  V0 = %f INITDUR = %f IHOLD = %f\n", V0, INITDUR, IHOLD)

  voltageClamp.rs = 1e9 // so the current it delivers during a run is miniscule
    // this is a "suspenders & belt" approach because Vsource[0].toff = 0
    // should prevent it from delivering nonzero current when t>0.
  *****/

  // restore simulation parameters
  dt = dtsav
  tstop = tstopsav
  t = 0

  // restore and re-init cvode if necessary
  if (temp!=0) {
    cvode.active(1)
    cvode.re_init()
  } else {
    fcurrent()
  }
  frecord_init()
  Log(LOG_TRACE, "Leaving method init()\n")
}

/******** End functions from custominit.hoc ***********/

/**************************************************
 * The following function is based on the contents of Christina's
 * simulateCurrentInj.hoc script
 */

/*******************
	simulateCurrentStep_withIHold

	Take the model with its currently defined parameters, and simulate a current step of
	specified size.  Uses the function eval_FRandCV(), found in rec_volt_justV.hoc

	Arguments:
	  $1	time at which the current step begins
		$2	time at which the current step ends
		$3	amount of current to inject, in nA (includes holding current)
		$4	1 or 0:  create binary file with (t,V) data?
		$s5	file name for output (see readNRNbin_Vonly.m to read this in MATLAB)
*******************/
proc simulateCurrentStep_withIHold() { local firingRate, startWindow, endWindow
  Log(LOG_TRACE, "Entering method simulateCurrentStep_withIHold()\n")

  ihold.del = 0
  ihold.dur = 1e9
  // Make sure that IHOLD has been computed.  Note that, in order for this to
  // work, ihold the IClamp must have already been set up to have an 'infinite'
  // duration, because it's used in the computation of IHOLD the value.  (The
  // amplitude of ihold the IClamp is then set to match IHOLD the value in the
  // code below.)
  init()

  ihold.amp = IHOLD
//  print "Done ihold = ", IHOLD

  stim.del = $1
  stim.dur = $2

  startWindow = $1 + 200  // Just for testing
  endWindow = $2
//  firingRate = eval_FRandCV($1, $2, $3 - ihold.amp, $s5, $4, 1)
  firingRate = eval_FRandCV($1, $2, $3, $s5, $4, 1)
  printf("Mean firing rate: %f\n", firingRate)
  Log(LOG_TRACE, "Leaving method simulateCurrentStep_withIHold()\n")
}

/******** End functions from simulateCurrentInj.hoc ***********/

/*
 * For the specified neuron, applies a stim current of -50pA to 50pA, at
 * 10pA intervals, in addition to a holding current that would by itself
 * maintain the neuron at a standardized potential specified by the global
 * V0.  At each stim current, the simulation progresses to steady state, which
 * is assumed to occur after 1000 ms, and then the membrane voltage and input
 * current are printed.  This produces a table of data that can be fit to
 * compute the input resistance of the cell.
 * (Procedure from Luebke and Rosene, J. Comp. Neurol. 2003)
 *
 * Arguments:
 */
proc printInputResistanceValuesForNeuron() { local test_amp
  stim.amp = 0

  // Set up the ihold current clamp so it can be used to compute the holding
  // current.
  ihold.del = 0
  ihold.dur = 1e9

  // Compute the holding current as a side effect.  The result is stored in the
  // global IHOLD
  init()
  ihold.amp = IHOLD

  stim.del = 100
  stim.dur = 1100

  tstop = 1100

  strdef test_amp_string
  // Neuron's IClamp is measured in nA, the test current is in pA.
  for (test_amp = -0.05; test_amp < 0.06; test_amp = test_amp + 0.01) {
    sprint(test_amp_string, "%.2f", test_amp)
    stim.amp = test_amp
    volt2txt(test_amp_string)
    printf("%s\t%f\n", test_amp_string, soma.v(0.5))
  }

}


proc sim_all_FR() { local fval, iinj, stimStop, cvval

  stim.del = 215
  stimStop = 2015
  stimStop = $1

  for( iinj = .13; iinj <= .33; iinj = iinj + .1 ) {
      fval = eval_FRandCV(215,stimStop,iinj,"", 0)
    if( fr.size > 1 ) { cvval = fr.stdev / fr.mean } else {cvval = 0 }
    printf("\tInject %g pA \t FR = %g Hz\n",iinj,fval) 
    //if( FRout.isopen() ) FRout.printf("%g\t%g\t%g\t",iinj,fval,cvval)
  }
}

proc sim_RN() {  local iinj, fval

  stim.del = 500
  stim.dur = 400
  for( iinj = -.04; iinj <= 0; iinj = iinj + .02 ) {
      fval = eval_FRandCV(500,900,iinj,"", 0)
    //printf("%g pA = %g Hz\n",1/g_pas,e_pas,gbar_ar, iinj,fval)
    //FRout.printf("%g\t%g\t",iinj,soma.v(.5))
  }
}

xopen("scaleRm_aug3f.hoc")

proc simulateNeuron() { local fval, iinj, RmScaled, CmScaled, nidx, custCM, task_idx

  strdef neuron_path, neuron_out
  neuron_path = $s1
  nidx = $2
  custCM = $3
  task_idx = $4

  printf("*****Entering simulateNeuron()\n\tNeuron path: %s\n", neuron_path)

  IHOLD = 0

  Log(LOG_INFO, "About to read cell\n")
  readcell(neuron_path, 3, RA, CM)  // 3: apical and basal both present
  Log(LOG_INFO, "Done reading cell\n")

  objref stim
  soma {
    Log(LOG_INFO, "Creating both stim and ihold\n")
    stim = new IClamp(0.5)
    ihold = new IClamp(0.5)
  }
  // Per Christina; this prevents interference during initialization.
  stim.del = 1e9


/** don't use the voltageClamp for Aniruddha's model? Use the shunt instead.  **/

  /*************
  // Voltage clamp, used to determine the holding current necessary to keep the
  // neuron at V0.
  soma {
    voltageClamp = new Vsource(0.5)
  }
  voltageClamp.rs = 1  // Internal resistance, megaohm.
  voltageClamp.toff = 0  // Time at which the voltageClamp ceases
  voltageClamp.amp = 0  // Target voltage for the clamp (the variable name is a
                        // bit of a misnomer.
  *************/

  /*********  	Add a shunt, as Aniruddha does (implemented before we started using the 
  **		Vsource or SEClamp mechanisms.  
  ********/


  /**************************************************
  * The following is taken from the global initialization of Christina's
  * aux_procs.hoc script.
  */


  add_yadav_axon()
  setup_SecLists()

  load_file("linear_conductances_traub.hoc")

  init_yadav_model()
  set_dataVec()

  forall cm = CM
  forall Ra = RA
  forall e_pas = E_PAS
  objref ptbVec // Perturbation vector?


  /******** End global initialization of aux_procs.hoc ***********/

  // Ensure a single-compartment soma
  soma {
    nseg = 1
    // Ensures that all neurons have the same size soma, so that we're comparing
    // dendritic effects alone.
    L = STD_SOMA_RADIUS * 2
    diam = STD_SOMA_RADIUS * 2
  }

  /**************************************************
  * The following is taken from Christina's main_PFCwSEClamp_forPCoskren.hoc
  * script, with bits and pieces extracted from or set based on the scripts
  * rigPFCmod.ses and vsrc.ses.
  */
  // Parameters tuned by Christina Weaver, Nov 2011, to fit representative physiology of one Jennie's young PFC neurons.
  RmScaled = RMVAL * scaleRm_vsAug3f(nidx)
  if( custCM == 0 ) {
    CmScaled = CM 
    printf("\tNo scaling of baseline CM 0.833\n")
  } else {
    CmScaled = 0.6* CM * scaleCm_vsAug3f(nidx)
    printf("\tCM scaled from baseline 0.833 to %.3f\n",CmScaled)
  }
  set_memres(RmScaled)
  adj_Cm(CmScaled)
  set_dend_Hratios(HSLP)

  load_file("Vkeep.ses")
  if( task_idx == 0 ) { 
      printf("\tSimulating injections below AP threshold\n")
      sim_RN()
  } else {
      printf("\tSimulating several injections above AP threshold\n")
      sim_all_FR(1215)
  }

}


steps_per_ms=20
dt=.05

// PICK UP HERE - RUN THE FI CURVES!

FRout = new File()
HVAL = .0001
E_PAS = -70	// Aniruddha's value
HSLP = 1
RMVAL = 17498
NAF_RED = 1
CAVAL = 5.41105e-05 
KC_RED = 1
PHI_SCL = 1
BETA_SCL = 1


proc run_YadavTraub() { 
    simulateNeuron(neuron_paths.o($1).s(),$1,$2,$3)
}