Rhesus Monkey Layer 3 Pyramidal Neurons: Young vs aged PFC (Coskren et al. 2015)

 Download zip file   Auto-launch 
Help downloading and running models
Accession:168858
Layer 3 (L3) pyramidal neurons in the lateral prefrontal cortex (LPFC) of rhesus monkeys exhibit dendritic regression, spine loss and increased action potential (AP) firing rates during normal aging. The relationship between these structural and functional alterations, if any, is unknown. Computational models using the digital reconstructions with Hodgkin-Huxley and AMPA channels allowed us to assess relationships between demonstrated age-related changes and to predict physiological changes that have not yet been tested empirically. Tuning passive parameters for each model predicted significantly higher membrane resistance (Rm) in aged versus young neurons. This Rm increase alone did not account for the empirically observed fI-curves, but coupling these Rm values with subtle differences in morphology and membrane capacitance Cm did. The predicted differences in passive parameters (or other parameters with similar effects) are mathematically plausible, but must be tested empirically.
Reference:
1 . Coskren PJ, Luebke JI, Kabaso D, Wearne SL, Yadav A, Rumbell T, Hof PR, Weaver CM (2015) Functional consequences of age-related morphologic changes to pyramidal neurons of the rhesus monkey prefrontal cortex. J Comput Neurosci 38:263-83 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Neuron or other electrically excitable cell;
Brain Region(s)/Organism:
Cell Type(s): Neocortex V1 L2/6 pyramidal intratelencephalic GLU cell;
Channel(s): I Na,t; I A; I K; I M; I h; I K,Ca; I Calcium; I_AHP;
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: NEURON;
Model Concept(s): Influence of Dendritic Geometry; Detailed Neuronal Models; Action Potentials; Aging/Alzheimer`s;
Implementer(s): Weaver, Christina [christina.weaver at fandm.edu];
Search NeuronDB for information about:  Neocortex V1 L2/6 pyramidal intratelencephalic GLU cell; I Na,t; I A; I K; I M; I h; I K,Ca; I Calcium; I_AHP;
objref voltage_vec, time_vec
// load_file("fixnseg.hoc")
// load_file("readcell_nomechanisms.hoc")

// *** 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 = 3

// 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())
  }
}


// *** Globals ***

// Variables that define the parameter space from Christina's V1 paper.  These can be modified
// to sample firing rates at different points
NA_CONDUCTANCE = ${gNa}
K_CONDUCTANCE = ${gKv}
stifleCharCount = printf("Running simulation for point in parameter space: %f, %f\n", NA_CONDUCTANCE, K_CONDUCTANCE)

// 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

// 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 = ${stdSomaRadius}  // Only used for electrotonic measurements
// CM = 1
// RA = 150
// celsius = 21

stifleCharCount = printf("STD_SOMA_RADIUS = %f\n", STD_SOMA_RADIUS)

// Global variables from Christina's add_axon.hoc script.
n_axon_seg  	= 5		/* # nodes in synthetic axon */
create iseg, hill, myelin[n_axon_seg], node[n_axon_seg]
objref axExcit, axSame

// 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 = ${V0}

stifleCharCount = printf("V0 = %f\n", V0)

// 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_resting_potentials
objref neuron_path_str
objref neuron_name_str

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

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

  connect_axon()
  insert_channels()
  active_set()
  passive_set()
  Axon = 1
  define_shape()

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

proc connect_axon() { local a,i		// cmw Aug '11:  deleted local var 'n'
  Log(LOG_TRACE, "Entering method connect_axon\n")

	create iseg
	create node[n_axon_seg]
	create hill
	create myelin[n_axon_seg]

	a = 0
	soma {
    for(x) a += area(x)
		equiv_diam = sqrt(a/(4*PI))
  }

	for i=0,n_axon_seg-1 {
    iseg       { L = 15  nseg = 5  diam = equiv_diam / 10 } /*Sloper&Powell 1982,Fig.71*/
    myelin[i]  { L = 100 nseg = 5  diam = iseg.diam  }
    node[i]    { L = 1.0 nseg = 1  diam = iseg.diam*.75 }
  }
	hill  { L = 10  nseg = 5  diam(0:1) = 4 * iseg.diam:iseg.diam }

	soma connect hill(0), 0.5
	hill connect iseg(0), 1
	iseg connect myelin[0](0), 1
	myelin[0] connect node[0](0), 1
	for i = 0,n_axon_seg-2  {
		node[i]	connect myelin[i+1](0),1
		myelin[i+1] connect node[i+1](0),1
	}

	axExcit = new SectionList()
	axSame  = new SectionList()
	hill axExcit.append()
	iseg axExcit.append()
	forsec "myelin" axSame.append()

	Axon = 1

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

proc insert_channels() {	/* insert channels and set reversal potentials */
  Log(LOG_TRACE, "Entering method insert_channels\n")

  forall {
    insert pas  /* generic conductance with reverse potential */
    insert na
    insert kv
  }

  // Vetter et al allowed the option to include Q, Ca, KCa, KM currents here.
  // cmw:  DELETED Aug '11

  forsec "myelin" uninsert kv /* no delayed rectifiers in myelin */

  /* set reversal potentials */
  forall e_pas = E_PAS
  forall if(ismembrane("k_ion"))  ek  	    = Ek
  forall if(ismembrane("na_ion")) ena	      = Ena
  forall if(ismembrane("na_ion"))	vshift_na = -10.5
  forall if(ismembrane("ca_ion")) {
    eca = 140
    ion_style("ca_ion", 0, 1, 0, 0, 0)
    vshift_ca = 0
  }
  Log(LOG_TRACE, "Leaving method insert_channels\n")
}

/******** End functions from add_axon.hoc ***********/

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

proc init_model() {
  Log(LOG_TRACE, "Entering method init_model\n")
	get_standard()		// initialize standard settings  (T, activespine ...)
	act0_set()			// just use the Vetter et al default
	active_set()
	passive_set()
  Log(LOG_TRACE, "Leaving method init_model\n")
}

/*
 * Establishes globals for the simulations.
 *   What does "substitutively" mean in this context?
 */
proc get_standard() {
  Log(LOG_TRACE, "Entering method get_standard\n")
  Iq_current = 0     // substitutively no Iq_current
  Ca_current = 0     // substitutively no Ca_current
  KCa_current = 0    // substitutively no KCa_current
  Km_current = 0     // substitutively no Km_current
  nonuniform_Rm = 0  // substitutively no nonuniform Rm
  electrotonicL = 0  // substitutively physical lengths
  activespine = 1    // model spines with active membrane
  currentdist = 1    // status for which distlist to use
  currentcell = 0    // number of currently active cell
  simMode = 0        // 0 IClamp 1 Vclamp simulation
  originx = 0.5      // where to add axon
  Ek = -90
  Ena = 60
  cells = 1
  swc = 0            // flag is turned on for Duke Southampton files
  synthetic = 0      // synthetic neurons need longer simulation time

  // stimulation (duration,delay,amplitude)
  St_del = 0
  St_amp = 1
  St_dur = 9999

  // simulation duration and dt.  Run for Sim_durI with Sim_dtI, then Sim_dur
  // with Sim_dt
  Sim_durI = 0
  Sim_dtI = 0.25
  Sim_dur = 15
  Sim_dt = 0.025

  // recording options (for display)
  Rec_del = Sim_durI
  Rec_dur = Sim_dur
  Rec_dt = Sim_dt

  // passive membrane properties, which aren't being changed in general
  Rm_node = 50
  Cm_myelin = 0.04

  G_ca = .3  /* 13.8.98 */  // What does this mean?
  G_kca = .1
  G_km = 3
  Rm_end = 0
  Rm_half = 0
  Rm_steep = 0
  G_qq = 0.02  // Iq settings
  Qq_end = 20
  Qq_steep = 439
  Qq_half   = 50

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

/* Zach Nature standard (kvz|naz) */
proc act0_set() {
  Log(LOG_TRACE, "Entering method act0_set\n")

  Rax		= RA		/* passive membrane properties */
	Rm         	= RM
	C_m        	= CM
	G_na		= NA_CONDUCTANCE  /* active membrane properties */
	G_kv		= K_CONDUCTANCE
	G_nanode	= 35000
	G_kvnode	= 500
  insert_channels()

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

proc active_set()  {	/* set active model parameters  */
  Log(LOG_TRACE, "Entering method active_set\n")
  g_na     = G_na
  g_kv     = G_kv
  g_nanode = G_nanode
  g_kvnode = G_kvnode
  g_qq		 = G_qq
  qq_end	 = Qq_end
  qq_steep = Qq_steep
  qq_half	 = Qq_half

  if (Iq_current)	forall {
    gbar_qq = qq_end + \
        (g_qq - qq_end) / (1 + exp((distance(0) - qq_half) / qq_steep))
  }

  forsec "soma" {
    gbar_na	= g_na
    gbar_kv = g_kv
  }

  forsec "myelin" gbar_na = g_na

  forsec "node" {
    gbar_na = g_nanode
    gbar_kv = g_kvnode
  }

  forsec "hill" {
    gbar_na = g_nanode
    gbar_kv = g_kvnode
  }

  forsec "iseg" {
    gbar_na = g_nanode
    gbar_kv = g_kvnode
  }

  // Vetter et al used different settings here for spines vs. dendrites
  forsec "dend" {
    gbar_kv = g_kv*dendscale
    gbar_na = g_na*dendscale
    if (Iq_current) gbar_qq *= dendscale
  }
  Log(LOG_TRACE, "Leaving method active_set\n")
}

proc passive_set() { /* set passive properties	*/
  Log(LOG_TRACE, "Entering method passive_set\n")
  /* $$1 	rm		*/
  /* $$2 	rax		*/
  /* $$3 	c_m 		*/
  /* $$4 	cm_myelin	*/
  /* $$5 	qq_end	 	*/
  /* $$6 	qq_steep 	*/
  /* $$7 	qq_half	 	*/

  n = numarg()
  if (n>0)  rm        = $$1      else rm     	= Rm
  if (n>1)  rax       = abs($$2) else rax    	= Rax
  if (n>2)  c_m       = $$3      else c_m    	= C_m
  if (n>3)  cm_myelin = $$4      else cm_myelin	= Cm_myelin
  if (n>4)  rm_node   = $$5      else rm_node	= Rm_node
  if (n>5)  rm_end    = $$6      else rm_end	= Rm_end
  if (n>6)  rm_steep  = $$7      else rm_steep	= Rm_steep
  if (n>7)  rm_half   = $$8      else rm_half	= Rm_half

  if (!nonuniform_Rm) rm_end = rm	 /* only nonuniform distr if nonuniform_Rm==1 */
  forall {
    if (rm - rm_end == 0) {
      g_pas = 1 / rm
    } else {
      g_pas = 1 / (rm + (rm - rm_end) / \
                        (1 + exp((distance(0) - rm_half) / rm_steep)))
    }
    cm    = c_m
    Ra    = rax
  }
  forsec "myelin"  	cm    = cm_myelin
  forsec "node"    	g_pas = 1 / rm_node
  forsec "iseg"    	g_pas = 1 / rm	/* make sure iseg & hill have uniform rm */
  forsec "hill"    	g_pas = 1 / rm
  Log(LOG_TRACE, "Leaving method passive_set\n")
}

/******** End functions from init_model.hoc ***********/

/**************************************************
 * 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 recording file name ->%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()
	printf("Ihold -> del: %f dur: %f, amp: %f\n", ihold.del, ihold.dur, ihold.amp)
	printf("Stim -> del: %f dur: %f, amp: %f\n", stim.del, stim.dur, stim.amp)
	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()
  }
  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
  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
 */

// $$1: Basic value for gbar_na
// $$2: Multiplicative factor to apply to gbar_na in the axon hillock and initial segment.
proc setGNaWithAxonHillockMultiplier() {
  Log(LOG_TRACE, "Entering method setGNaWithAxonHillockMultiplier\n")
  soma gbar_na = $$1
  forsec "dend" gbar_na = $$1
  forsec axSame gbar_na = $$1
  forsec axExcit gbar_na = $$1 * $$2
  Log(LOG_TRACE, "Leaving method setGNaWithAxonHillockMultiplier\n")
}

// $$1: Basic value for gbar_kv
// $$2: Multiplicative factor to apply to gbar_kv in the axon hillock and initial segment.
proc setKvWithAxonHillockMultiplier() {
  Log(LOG_TRACE, "Entering method setKvWithAxonHillockMultiplier\n")
	forall {
		if(ismembrane("kv")) gbar_kv = $$1
	}
	forsec axExcit gbar_kv = $$1 * $$2
  Log(LOG_TRACE, "Leaving method setKvWithAxonHillockMultiplier\n")
}

proc set_gpas() {
  Log(LOG_TRACE, "Entering method set_gpas\n")
  forall {
    ifsec "node" continue
    g_pas = $$1
	}
  Log(LOG_TRACE, "Entering method set_gpas\n")
}

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

/******** End functions from aux_procs.hoc ***********/

/**************************************************
 * 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

	The amperage specified in parameter 3 is the *total* current that should be
	injected into the neuron, including both the holding current and the test
	current.  Thus, if argument 3 is 130 nA, and the holding current is 3 nA, then
	the "ihold" IClamp object will have an amp of 3, and the "stim" IClamp object
	will have an amp of 127.  Thus all neurons receive the same total current
	injection, even if they require different holding currents to bring their
	resting potential to a standard value.

	Arguments:
	  $$1	time at which the current step begins
		$$2	time at which the current step ends
		$$3	total amount of current to inject, in nA
		$$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 prepNeuronForSimulation() {
  strdef neuron_path
  neuron_path = $$s1

  printf("Neuron path: %s\n", neuron_path)

  IHOLD = 0

  Log(LOG_INFO, "About to read cell\n")
  readcell(neuron_path, 3)  // 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

  // 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.

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

  init_model()
  add_axon()
  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

    printf("Using soma radius of %f\n", STD_SOMA_RADIUS)
  }

  /**************************************************
  * 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.
  set_gpas(G_PAS)
  setGNaWithAxonHillockMultiplier(${gNa},1e3)
  setKvWithAxonHillockMultiplier(${gKv},16.667)
  set_epas(E_PAS)

}

prepNeuronForSimulation("${neuronPath}")
//  simulateCurrentStep_withIHold(100, 1100, 0, 1, "tmp")
// printInputResistanceValuesForNeuron(neuron_path)

// quit()

Loading data, please wait...