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