original.hoc

//Template file for all cell types used in the network described in Zhang TC et al., J. Neurophysiology, 2014.  Relevant geometric and cellular parameters are listed within each "template."  Units for all values used listed at bottom (scroll down) of .hoc file.

load_file("nrngui.hoc")

begintemplate MelnickSG //This is "SG Cell" and "SGSCS Cell"

//From Melnick et al. J. Physiol (2004) Modified by Tianhe C. Zhang (tz5@duke.edu).

ndend=1

public soma, hillock, axon, dend
public synlist, x, y, z, position, connect2target
objref synlist, syn_
create soma, hillock, axon, dend


proc init(){
access soma
{diam=10 L=10 nseg=10} // area 100 um2 means mA/cm2 identical to nA 
insert B_Na
insert B_A
insert B_DR
insert KDR
insert KDRI
insert pas
{gnabar_B_Na = 0.008 ena = 60 g_pas = 1.1e-05 e_pas = -70 gkbar_KDRI = 0.0043 ek = -84}

access hillock
{L=30 nseg=30 dsoma=1 daxon=0.5 diam(0:1)=dsoma:daxon} 
insert B_Na
insert B_A
insert B_DR
insert KDR
insert KDRI
insert pas
{gnabar_B_Na = 3.45 ena = 60 g_pas = 1.1e-05 e_pas = -70 gkbar_KDRI = 0.076 ek = -84} // Tonic firing with maintained I-F characteristics.

access axon //Stump to measure AP propagation.  No other purpose.
{diam=0.001 L=0.001 nseg=50}
insert B_Na
insert B_A
insert B_DR
insert KDR
insert KDRI
insert pas
{gnar_B_Na = 0 ena = 60 g_pas = 0 e_pas = -70 gkbar_KDRI = 0 ek = -84}

access dend
{nseg=50 diam=1.4 L = 1371}
insert SS
insert B_DR
insert KDR
insert KDRI	
insert pas
{gnabar_SS = 0 ena = 60 g_pas = 1.1e-05 e_pas = -70 gkbar_KDRI = 0.034 ek = -84}

soma connect hillock(0),1
hillock connect axon(0),1
soma connect dend(0),0

forall Ra = 80 

synlist = new List()

synapses()

x = y = z = 0

}

objref SynapseID, SynapseFile, NumVec
strdef SynapseIDFile

proc synapses(){
//Uploads list of files that contains the number of each type of synapse (see below) corresponding to each neuron of this type in the system.

SynapseFile = new File()
SynapseID = new File()
SynapseID.ropen("SG_Cell_SynapseID.dat")
SynapseID.scanstr(SynapseIDFile)
SynapseID.close()
NumVec = new Vector()

SynapseFile.ropen(SynapseIDFile)
NumVec.scantil(SynapseFile, -1e15)
SynapseFile.close()

//Assigns synapses based on synapse counts obtained from "SynapseFile".  Unfortunately, synapse composition can't be dynamically defined using HOC, so one must redo the list below if synapse counts need to be changed.

for num = 0, NumVec.x[0]-1{
dend syn_ = new AMPA_DynSyn(0.5) syn_.tau_rise = 0.1 syn_.tau_decay = 5 synlist.append(syn_)
}
for num = 0, NumVec.x[1]-1{
dend syn_ = new AMPA_DynSyn(0.5) syn_.tau_rise = 0.1 syn_.tau_decay = 5 synlist.append(syn_)
}
}

proc position() { local i  // Optional feature.  This allows coupling of neuron to extracellular space e.g. for simulations of extracellular stimulation.
  soma for i = 0, n3d()-1 {
    pt3dchange(i, $1-x+x3d(i), $2-y+y3d(i), $3-z+z3d(i), diam3d(i))
  }
  x = $1  y = $2  z = $3
}

obfunc connect2target() { localobj nc //$o1 target point process, optional $o2 returned NetCon.  "Black box."
  soma nc = new NetCon(&v(1), $o1)
  nc.threshold = -30
  if (numarg() == 2) { $o2 = nc } // for backward compatibility
  return nc
}

endtemplate MelnickSG


//*NEW TEMPLATE*//
//*NEW TEMPLATE*//
//*NEW TEMPLATE*//


begintemplate AguiarIN// This is "EX Cell".

//Created by Paulo Aguiar [pauloaguiar@fc.up.pt].  Modified by Tianhe C. Zhang (tz5@duke.edu)

public soma, dend, hillock, axon, synlist, x, y, z, position, connect2target
create soma, dend, hillock, axon
objref synlist, syn_

proc init() {
x = y = z = 0
create soma // Soma
 
     soma {    
	nseg = 3  
	L = 20.0
	diam = 20.0
	
	//HH channels: iK
	insert HH2 {
	    gnabar_HH2 = 0.00  // Na not in soma.
	    gkbar_HH2 = 0.0043*1/4 // Scaled to safronov et al. (1997) and Wolff et al.(1998) ratios.
	    vtraub_HH2 = -55.0
	}
	
	//intracellular Ca dynamics
	insert CaIntraCellDyn {
	    depth_CaIntraCellDyn = 0.1
	    cai_tau_CaIntraCellDyn = 1.0
	    cai_inf_CaIntraCellDyn = 50.0e-6
	}
	
	//potassium current dependent on intracellular calcium concentration 
	insert iKCa {
	    gbar_iKCa = 0.002 //0.002
	}
	
	insert pas // Passive parameters.
	g_pas = 4.2e-5
	e_pas = -65.0
	ek = -70.0
	Ra = 150.0
    }
    
create dend // Single dendrite. 
 
    dend {    
	nseg = 5    
	L = 400.0  
	diam = 3.0

	//HH channels: iNat and iK
	insert HH2{
	    gnabar_HH2 = 0.000
	    gkbar_HH2 = 0.036 // *5/6
	}
	
	//intracellular Ca dynamics
	insert CaIntraCellDyn {
	    depth_CaIntraCellDyn = 0.1
	    cai_tau_CaIntraCellDyn = 2.0
	    cai_inf_CaIntraCellDyn = 50.0e-6
	}

	//potassium current dependent on intracellular calcium concentration 
	insert iKCa {
	    gbar_iKCa = 0.002 //0.002
	}	
	
	insert pas
	g_pas = 4.2e-5
	e_pas = -65.0
	ek = -70.0
	Ra = 150.0
    }
    
    
create hillock

    hillock {   
	nseg = 3  
	L = 9.0
	diam(0:1) = 2.0:1.0
	
	//HH channels: iNa,t and iK
	insert HH2 {
	    gnabar_HH2 = 3.45
	    gkbar_HH2 = 0.076
	    vtraub_HH2 = -55.0
	}
	
	insert pas
	g_pas = 4.2e-5
	e_pas = -65.0
	Ra = 150.0
    }
    
    create axon
    axon {    
	nseg = 5
	L = 1000.0
	diam = 1.0
	
	//HH channels: iNa,t and iK
	insert HH2 {
	    gnabar_HH2 = 0.0
	    gkbar_HH2 = 0.00	//0.06
	    vtraub_HH2 = -55
	}
	
	insert pas
	g_pas = 4.2e-5
	e_pas = -65.0
	Ra = 150.0
    }
    
    //CONNECTIONS between cell elements.
    soma connect hillock(0),1
    hillock connect axon(0),1
    soma connect dend(0),0 

    synlist = new List()
    synapses()
}

proc position() { local i // Optional feature.  This allows coupling of neuron to extracellular space e.g. for simulations of extracellular stimulation.
  soma for i = 0, n3d()-1 {
    pt3dchange(i, $1-x+x3d(i), $2-y+y3d(i), $3-z+z3d(i), diam3d(i))
  }
  x = $1  y = $2  z = $3
}

obfunc connect2target() { localobj nc //$o1 target point process, optional $o2 returned NetCon
  soma nc = new NetCon(&v(1), $o1)
  nc.threshold = -30
  if (numarg() == 2) { $o2 = nc } // for backward compatibility
  return nc
}


objref SynapseID, SynapseFile, NumVec
strdef SynapseIDFile

proc synapses(){
//Uploads list of files that contains the number of each type of synapse (see below) corresponding to each neuron of this type in the system.

SynapseFile = new File()
SynapseID = new File()
SynapseID.ropen("EX_Cell_SynapseID.dat")
SynapseID.scanstr(SynapseIDFile)
SynapseID.close()
NumVec = new Vector()

SynapseFile.ropen(SynapseIDFile)
NumVec.scantil(SynapseFile, -1e15)
SynapseFile.close()

//Assigns synapses based on synapse counts obtained from "SynapseFile".  Unfortunately, synapse composition can't be dynamically defined using HOC, so one must redo the list below if synapse counts need to be changed.

for num = 0, 29{
dend syn_ = new AMPA_DynSyn(0.5) syn_.tau_rise = 0.1 syn_.tau_decay = 5 synlist.append(syn_)
}
for num = 0, 29{
dend syn_ = new NMDA_DynSyn(0.5) syn_.tau_rise = 2 syn_.tau_decay = 100 synlist.append(syn_)
}
for num = 0, 29{
dend syn_ = new NK1_DynSyn(0.5) syn_.tau_rise = 100 syn_.tau_decay = 3000 synlist.append(syn_)
}
for num = 0, 1{
dend syn_ = new GABAa_DynSyn(0.5) syn_.tau_rise = 0.1 syn_.tau_decay = 20 syn_.e = -70 synlist.append(syn_)
}
for num = 0, 1{
dend syn_ = new GABAa_DynSyn(0.5) syn_.tau_rise = 0.1 syn_.tau_decay = 20 syn_.e = -70 synlist.append(syn_)
}

}

endtemplate AguiarIN


//*NEW TEMPLATE*//
//*NEW TEMPLATE*//
//*NEW TEMPLATE*//

//Created by Paulo Aguiar [pauloaguiar@fc.up.pt].  Modified by Tianhe C. Zhang (tz5@duke.edu)

// CREATE WDR NEURON

begintemplate AguiarWDR // This is "T_Cell"

public soma, dend, hillock, axon, x, y, z, position, connect2target
public synlist
objref syn_, synlist
create soma, dend, hillock, axon

proc init() {
 
x = y = z = 0
cascale = 1 // Scale calcium currents responsible for wind-up.  Set to 0 to eliminate calcium currents without having to change template codes below.  Do not set to a negative value unless a very good reason is found to do so.

create soma    

    soma {    
	nseg = 3  
	L = 20.0
	diam = 20.0
	
	//HH channels: iNat and iK
	insert HH2 { //No AP Initiation in the soma?
		gnabar_HH2 = 0.000 //"Default" 0.08.  Melnick 0.008.
		gkbar_HH2 = 0.0043*6/24  //0.3e-3 in P&D. "/4" is there in parallel with Safronov.
		vtraub_HH2 = -55.0
	}
	
	//intracellular Ca dynamics
	insert CaIntraCellDyn {
	    depth_CaIntraCellDyn = 0.1
	    cai_tau_CaIntraCellDyn = 1.0
	    cai_inf_CaIntraCellDyn = 50.0e-6
	}

	//high-voltage activated long-lasting calcium current, L-type
	insert iCaL {
	    pcabar_iCaL = 0.0001 * cascale//0.0001- IMPORTANT: this current drives the (activity control) somatic iKCa current
	}      
 
	//non-specific current dependent on intracellular calcium concentration 
	insert iCaAN {
	    gbar_iCaAN = 0.0000 //0.00000
	}    

	//potassium current dependent on intracellular calcium concentration 
	insert iKCa {
	    gbar_iKCa = 0.0001 * cascale //0.0001
	}

	//sodium persistent current
	insert iNaP {
	    gnabar_iNaP = 0.0001 * cascale//0.0001
	}    

	insert pas
	g_pas = 4.2e-5
	e_pas = -65.0
	ek = -70.0
	Ra = 150.0
    }
    
create dend

    dend {    
	nseg = 5    
	L = 350.0     
	diam = 2.5
	
	//intracellular Ca dynamics
	insert CaIntraCellDyn {
	    depth_CaIntraCellDyn = 0.1
	    cai_tau_CaIntraCellDyn = 2.0
	    cai_inf_CaIntraCellDyn = 50.0e-6
	}

	//high-voltage activated long-lasting calcium current, L-type
	insert iCaL {
	    pcabar_iCaL = 3e-5 * cascale //3.0e-5 IMPORTANT: this current is important for activity control (drives the iKCa current).  cascale currently set to 1.
	}             
    
	//non-specific current dependent on intracellular calcium concentration 
	insert iCaAN {
	    gbar_iCaAN = 0.00007 * cascale * 1.3 //0.00007; This is a sensible parameter
	    //higher values of gbar_iCaAN produce graphs similar to Silviloti et al 93.  Need 1.3 scalar to match up to Herrero (2000).
	}        

	//potassium current dependent on intracellular calcium concentration 
	insert iKCa {
	    gbar_iKCa = 0.001 * cascale//0.001; higher values place "holes" in the scatter plot, resulting from the cai bump after synaptic activation);
		//naturally, lower values will lead to increased firing
	}
	
	insert HH2{
	    gnabar_HH2 = 0.000
	    gkbar_HH2 = 0.036 // *5/6
	}

	insert pas
	g_pas = 4.2e-5
	e_pas = -65.0
	ek = -70.0
	Ra = 150.0
    }
    
create hillock

hillock {   
	nseg = 3  
	L = 9 //Note that default is 3.
	diam(0:1) = 2.0:1.0
	
	//HH channels: iNa,t and iK
	insert HH2 {
	    gnabar_HH2 = 3.45 //Remember Melnick is 3.45.  "/4" is for surface area compensation.
	    gkbar_HH2 = 0.076 // Default 0.04.  POss. need to scale by 0.25.  Check I-F with normal K.
	    vtraub_HH2 = -55.0
	}

	insert B_A{
	   gkbar_B_A = 0.0 // A-current for delayed firing.
	}

	insert pas
	g_pas = 4.2e-5
	e_pas = -65.0
	Ra = 150.0
    }
    
create axon

    axon {    // Axon is a single node "stump" used only to check action potential propagation.
	nseg = 5
	L = 1000
	diam = 1.0
	
	//HH channels: iNa,t and iK
	insert HH2 {
	    gnabar_HH2 = 0 //0.1
	    gkbar_HH2 = 0	//0.04
	    vtraub_HH2 = -55
	}

	insert pas // Axon is a single node "stump."
	g_pas = 0
	e_pas = -65.0	
	Ra = 150.0
    }
    
//CONNECTIONS between neural elements.
soma connect hillock(0),1
hillock connect axon(0),1
soma connect dend(0),0 

synlist = new List()
synapses()

}

objref SynapseID, SynapseFile, NumVec
strdef SynapseIDFile

proc synapses(){

//Uploads list of files that contains the number of each type of synapse (see below) corresponding to each neuron of this type in the system.

SynapseFile = new File()
SynapseID = new File()
SynapseID.ropen("T_Cell_SynapseID.dat")
SynapseID.scanstr(SynapseIDFile)
SynapseID.close()
NumVec = new Vector()

//Assigns synapses based on synapse counts obtained from "SynapseFile".  Unfortunately, synapse composition can't be dynamically defined using HOC, so one must redo the list below if synapse counts need to be changed.

SynapseFile.ropen(SynapseIDFile)
NumVec.scantil(SynapseFile, -1e15)
SynapseFile.close()

for num = 0, 14{
dend syn_ = new AMPA_DynSyn(0.5) syn_.tau_rise = 0.1 syn_.tau_decay = 5 synlist.append(syn_)
}
for num = 0, 14{
dend syn_ = new NMDA_DynSyn(0.5) syn_.tau_rise = 2 syn_.tau_decay = 100 synlist.append(syn_)
}
for num = 0, 14{
dend syn_ = new AMPA_DynSyn(0.5) syn_.tau_rise = 0.1 syn_.tau_decay = 5 synlist.append(syn_)
}
for num = 0, 14{
dend syn_ = new NMDA_DynSyn(0.5) syn_.tau_rise = 2 syn_.tau_decay = 100 synlist.append(syn_)
}
for num = 0, 29{
dend syn_ = new NK1_DynSyn(0.5) syn_.tau_rise = 200 syn_.tau_decay = 3000 synlist.append(syn_)
}
for num = 0, 29{
dend syn_ = new AMPA_DynSyn(0.5) syn_.tau_rise = 0.1 syn_.tau_decay = 5 synlist.append(syn_)
}
for num = 0, 29{
dend syn_ = new NMDA_DynSyn(0.5) syn_.tau_rise = 2 syn_.tau_decay = 100 synlist.append(syn_)
}
for num = 0, 0{
dend syn_ = new Glycine_DynSyn(0.5) syn_.tau_rise = 0.1 syn_.tau_decay = 10 synlist.append(syn_)
}
for num = 0, 0{
dend syn_ = new GABAa_DynSyn(0.5) syn_.tau_rise = 0.1 syn_.tau_decay = 20 syn_.e = -70 synlist.append(syn_)
}
for num = 0, 0{
dend syn_ = new GABAa_DynSyn(0.5) syn_.tau_rise = 0.1 syn_.tau_decay = 20 syn_.e = -70 synlist.append(syn_)
}
for num = 0, 0{
dend syn_ = new GABAa_DynSyn(0.5) syn_.tau_rise = 0.1 syn_.tau_decay = 20 syn_.e = -70 synlist.append(syn_)
}
}

proc position() { local i 
  soma for i = 0, n3d()-1 {
    pt3dchange(i, $1-x+x3d(i), $2-y+y3d(i), $3-z+z3d(i), diam3d(i))
  }
  x = $1  y = $2  z = $3
}

obfunc connect2target() { localobj nc //$o1 target point process, optional $o2 returned NetCon
  soma nc = new NetCon(&v(1), $o1)
  nc.threshold = -30
  if (numarg() == 2) { $o2 = nc } // for backward compatibility
  return nc
}

endtemplate AguiarWDR

begintemplate S_NetStim 

//DUMMY NetStims; These are only included because Syn objects must be connected to NetCon object as noted in General Readme file.

public pp, connect2target, x, y, z, position, is_art
objref pp
proc init() {
  pp = new NetStim()
    pp.interval = 0
    pp.number = 0
    pp.start = 0
}

func is_art() { return 1 }
obfunc connect2target() { localobj nc
  nc = new NetCon(pp, $o1)
  if (numarg() == 2) { $o2 = nc }
  return nc
}

proc position(){x=$1  y=$2  z=$3}
endtemplate S_NetStim


//Network specification interface.  This is called in "MakeNetwork.hoc."  Black box.
objref cells, nclist, netcon
{cells = new List()  nclist = new List()}

func cell_append() {cells.append($o1)  $o1.position($2,$3,$4)
	return cells.count - 1
}

func nc_append() {//srcindex, tarcelindex, synindex
  if ($3 >= 0) {
    netcon = cells.object($1).connect2target(cells.object($2).synlist.object($3))
    netcon.weight = $4   netcon.delay = $5
  }else{
    netcon = cells.object($1).connect2target(cells.object($2).pp)
    netcon.weight = $4   netcon.delay = $5
  }
  nclist.append(netcon)
  return nclist.count - 1
}



//************************************************************************************
//UNITS	    
//Category          Variable        [Units]  (notes)

//Time                  t	    [ms]
//Voltage               v	    [mV]
//Current               i	    [mA/cm2] (distributed)	[nA] (point process)
//Concentration		ko, ki, etc.	[mM]
//Specific capacitance	cm          [uf/cm2]
//Length                diam, L     [um]
//Conductance		g	    [S/cm2]  (distributed)	[uS] (point process)
//Cytoplasmic resistivity  Ra	    [ohm cm]
//Resistance	        Ri	    [10E6 ohm]

//Modified by Tianhe Zhang from NEURON Book Chapter 11 (Carnevale and Hines) to import necessary values from connectivity files.
//NOTE: For the "import spike times" simulation, spike times will be imported IN THE SIMULATION RUNNER ROUTINE.

//This script is generally a black box i.e. nothing needs to be changed here when only altering the network/cell properties.  However, if a new parameter needs to be incorporated
//(e.g. adding an "ERevVector.txt" import to run KCC2 simulations; see "importcritvalues()"), this should be altered to reflect changes.

////////// import critical values ///////////////// Alteration from NEURON code earlier.  This will allow for parameters to be imported from MATLAB.

objref InCells, InX, InY, InZ, InTypes, InFrom, InTo, InSynapse, InWeights, InDelays, InThresh
objref CellsVec, XVec, YVec, ZVec, CellTypeVec, FromVec, ToVec,  SynapseVec, WeightsVec, DelaysVec, ThreshVec


proc importcritvalues() { //These values (which define the cells and connectivities) of the neural network MUST be imported.

	//Cells//
	InCells = new File()
	CellsVec = new Vector()
	InCells.ropen("CellVector.txt")
	CellsVec.scanf(InCells)
	InCells.close()

	//Cell Coordinates//
	InX = new File()
	XVec = new Vector()
	InX.ropen("XVector.txt")
	XVec.scanf(InX)
	InX.close()

	InY = new File()
	YVec = new Vector()
	InY.ropen("YVector.txt")
	YVec.scanf(InY)
	InY.close()

	InZ = new File()
	ZVec = new Vector()
	InZ.ropen("ZVector.txt")
	ZVec.scanf(InZ)
	InZ.close()

	//Cell Types//  CellTypeVec contains a list of indices corresponding to the indices of CellsVec where cell types switch.

	InTypes = new File()
	CellTypeVec = new Vector()
	InTypes.ropen("CellTypeVector.txt")
	CellTypeVec.scanf(InTypes)
	InTypes.close()


	//Connections//
	InFrom = new File()
	FromVec = new Vector()
	InFrom.ropen("FromVector.txt")
	FromVec.scanf(InFrom)
	InFrom.close()

	InTo = new File()
	ToVec = new Vector()
	InTo.ropen("ToVector.txt")
	ToVec.scanf(InTo)
	InTo.close()

	InSynapse = new File()
	SynapseVec = new Vector()
	InSynapse.ropen("SynapseVector.txt")
	SynapseVec.scanf(InSynapse)
	InSynapse.close()

	//
}

importcritvalues()

//Import data regarding cell parameters (e.g. tau, intervals).  Not importing a value means default in NEURON is used, so optional (if value doesn't matter).

//Import data regarding connection parameters (weights, delays, synaptic thresholds).  Not importing a value means default in NEURON is used.

proc importWeights(){
	InWeights = new File()
	WeightsVec = new Vector()
	InWeights.ropen("WeightVector.txt")
	WeightsVec.scanf(InWeights)
	InWeights.close()
}

proc importDelays(){
	InDelays = new File()
	DelaysVec = new Vector()
	InDelays.ropen("DelayVector.txt")
	DelaysVec.scanf(InDelays)
	InDelays.close()
}

proc importThresh(){
	InThresh = new File()
	ThreshVec = new Vector()
	InThresh.ropen("ThresholdVector.txt")
	ThreshVec.scanf(InThresh)
	InThresh.close()
}

//proc importERev(){  

//Reversal potentials for synapses.  Commented out for Wind-Up simulations, but essential to reproduce KCC2 results.  Make sure synapses that use same ion
//have same reversal potential (KCC2: rest potentials unaffected).

//	InERev = new File()
//	ERevVec = new Vector()
//	InERev.ropen("ERevVector.txt")
//	ERevVec.scanf(InERev)
//	InERev.close()
//}

importWeights()
importDelays()
importThresh()

//Greater than equal to.

func ge() {
  if ($1<$2) {
    $1=$2
  }
  return $1
}

////////// create a network [NEURON generated] //////////
// argument is desired number of cells

objref CellTempFile
strdef cmd

proc createnet() { local p, q //Input 1: # of Cells.  Input 2: # of connections (not categorized yet.  NOTE THAT USER MUST MANUALLY CHANGE "cell_append" ACCORDING TO HOW MANY TYPES OF CELLS THERE ARE)
  $1 = ge($1,2) // force net to have at least two cells
  ncell = $1
  // so we can make a new net without having to exit and restart
  nclist.remove_all()
  cells.remove_all()
  CellTempFile = new File()
  CellTempFile.ropen("CellTempList.dat")

  for p = 1, CellTypeVec.size()-1{
	  CellTempFile.gets(cmd) //Scan in relevant cell template command
	  for q = CellTypeVec.x[p-1], CellTypeVec.x[p]-1{
		j = q
		execute1(cmd) //Run object generation command once for each relevant cell.  See "CellTempList.dat" for exact commands (list of "VARIABLE = new CELLCLASS").
	  }
  }
  CellTempFile.close()

  for p=0, $2-1 {  //Connect cells that were created above based on the layout described in the Connections Excel file.
    nc_append(FromVec.x[p], ToVec.x[p], SynapseVec.x[p], WeightsVec.x[p], DelaysVec.x[p]) //CELL COUNT STARTS AT ZERO.
	print p
  }
  objref netcon  // leave no loose ends (see nc_append())
  strdef cmd     // same as above, but for strings.
}

////////// specify parameters for each individual connection in network//////////
// call this settau() to avoid conflict with scalar tau

proc settau() { local i
	temp = $2
    	if (temp>0) {
    	cells.object($1).pp.tau = $2
   	 }
}

// reworked for individual cells/connections.
proc interval() { local i
	temp = $2
    	if (temp>0) {
   	 cells.object($1).pp.invl = $2
    	}
}

proc weight() { local i
    nclist.object($1).weight = $2
}

proc delay() { local i
  $2 = ge($2,0)  // min del is 0 ms
  del = $2
    nclist.object($1).delay = $2
}

proc thresh() { local i
    nclist.object($1).threshold = $2
}


////////// actually execute the above.  Easier if encapsulated into a single process //////////
createnet(CellsVec.size(), FromVec.size())

proc ConnParamSetter(){ local k
	for k=0, FromVec.size()-1{
		weight(k, WeightsVec.x[k])
		delay(k, DelaysVec.x[k])
		thresh(k, ThreshVec.x[k])
	}
}
ConnParamSetter()


load_file("nrngui.hoc")
//load_file("MakeNetMATLAB.hoc")

objref netcon, vec, spikes, nil, graster

objref SpikeStatsFile, numspikes


proc TempInit(){ //Temperature and time step.
	celsius = 36 //According to Aguiar
	dt = 0.0125 // Based on 1/8 of fastest tau in system "rule of thumb."
}

TempInit()


proc SynapseInit(){ //Initializes synapses by specifying the number of synapses.  Also sets tstop.
	numspikes = new Vector()
	SpikeStatsFile = new File()
	SpikeStatsFile.ropen("SpikeStatsVector.txt")
	numsources = SpikeStatsFile.scanvar()
	numspikes.scanf(SpikeStatsFile, numsources)
	tstop = SpikeStatsFile.scanvar()
	SpikeStatsFile.close()

}

SynapseInit()

objref SpikeVectorFile, SpikeVector[numspikes.size()]

proc SynapseLoad(){ local k //Loads synapse event i.e. input spike times into vectors, which will be themselves inserted into the appropriate synapse using the netcon.event command.

	SpikeVectorFile = new File()
	SpikeVectorFile.ropen("SpikeTimesVector.txt")
	for k=0, numspikes.size()-1{
		SpikeVector[k] = new Vector(numspikes.x[k], 0)
		SpikeVector[k].scantil(SpikeVectorFile, -1e15)
	}
	SpikeVectorFile.close()
}

SynapseLoad()

objref fih

//Required according to NEURON Book Chapter 11 to initialize cells for "stdinit" call and to load spike times.  Black Box.

fih = new FInitializeHandler("loadqueue()")

proc loadqueue() { local j
	for k = 0, numspikes.size()-1 {
		for j=0, numspikes.x[k]-1{
			nclist.object(k).event(SpikeVector[k].x[j])
		}
	}
}

proc init() { //Taken from Ch. 11 of Neuron book.  Initializes the solver.
	finitialize(v_init)
	if (cvode.active()) {
		cvode.re_init()
	} else {
		fcurrent()
	} frecord_init()
}

//Actually run simulations
proc run() {
	stdinit()
	continuerun(tstop)
}

objref OutputVectorsSG, OutputVectorsSGM, OutputVectorsSGH, OutputTimesSG, netconSG, nil
objref OutputVectorsZ2, OutputTimesZ2, netconZ2
objref OutputVectorsT, OutputVectorsTM, OutputVectorsTH, OutputTimesT, netconT
objref OutputVectorsEX, OutputVectorsEXM, OutputVectorsEXH, OutputTimesEX, netconEX


proc OutputData(){ local i //Outputs relevant values into variables specified 2 lines above.  Comment/uncomment what is needed; currently includes vm, spike times, m, n, h for each class of neuron (SG, EX, T) (only Vm and spike times are currently written to a file).

	OutputVectorsSG = new Vector() //needs to be modified to accomodate multiple cells.
	//OutputVectorsSGM = new Vector()
	//OutputVectorsSGH = new Vector()
	OutputVectorsSG.record(&MelnickSG[0].soma.v(0.5))
	//OutputVectorsSGM.record(&MelnickSG[0].soma.m_hh(0.5))
	//OutputVectorsSGH.record(&MelnickSG[0].soma.h_hh(0.5))
	OutputTimesSG = new Vector()
	MelnickSG[0].soma netconSG = new NetCon(&v(0.5), nil)
	netconSG.threshold = -30 //-30 mV AP threshold.
	netconSG.record(OutputTimesSG)
	objref netconSG

	OutputVectorsZ2 = new Vector() //needs to be modified to accomodate multiple cells.
	OutputVectorsZ2.record(&MelnickSG[1].soma.v(0.5))
	OutputTimesZ2 = new Vector()
	MelnickSG[1].soma netconZ2 = new NetCon(&v(0.5), nil)
	netconZ2.threshold = -30 //-30 mV AP threshold.
	netconZ2.record(OutputTimesZ2)
	objref netconZ2	

	OutputVectorsT = new Vector() //needs to be modified to accomodate multiple cells.
	//OutputVectorsTM = new Vector()
	//OutputVectorsTH = new Vector()
	OutputVectorsT.record(&AguiarWDR[0].soma.v(0.5))
	//OutputVectorsTM.record(&AguiarWDR[0].soma.m_hh(0.5))
	//OutputVectorsTH.record(&AguiarWDR[0].soma.h_hh(0.5))
	OutputTimesT = new Vector()
	AguiarWDR[0].soma netconT = new NetCon(&v(0.5), nil)
	netconT.threshold = -30 //-30 mv AP threshold.
	netconT.record(OutputTimesT)
	objref netconT

	OutputVectorsEX = new Vector() //needs to be modified to accomodate multiple cells.
	OutputVectorsEX.record(&AguiarIN[0].soma.v(0.5))
	OutputTimesEX = new Vector()
	AguiarIN[0].soma netconEX = new NetCon(&v(0.5), nil)
	netconEX.threshold = -30 //-30 mv AP threshold.
	netconEX.record(OutputTimesEX)
	objref netconEX

init()
run()
}

proc Step(){ //Use something like this to probe internal variables, except replace "print" command with record command of some sort.  Currently not called, but can be useful.
	index = 1
	finitialize(v_init)
	while (t<tstop) {
	print MelnickSG[1].soma.v(0.5)
	fadvance()
	}
}


//Outputs.  Saves results of "OutputData()" into relevant .dat files. .dat files are saved in 'double' (MATLAB: floating point, 64 bit) binary format.

objref DestFileSG, TimeFileSG, ReadFileSG, ReadTimeFileSG, TempVectorSG //ADDITION: Time File for recording of spike times out of different SG Cells.
strdef SaveToMeSG

proc saveSG(){ local k
	DestFileSG = new File()
	ReadFileSG = new File()
	TimeFileSG = new File() // File for AP Time vector storage.
	TempVectorSG = new Vector()
	ReadFileSG.ropen("SG_Cell_filenames.dat") //MAKE SURE THIS NAME MATCHES YOUR LIST OF OUTPUT FILES.	
	
	//for k=0, 2 { Save multiple neurons. OPTIONAL.  Will need to add bracketed index entry to TempVector (e.g. "TempVector[k]").  
		TempVectorSG = OutputVectorsSG
		ReadFileSG.scanstr(SaveToMeSG)
		DestFileSG.wopen(SaveToMeSG)
		TempVectorSG.fwrite(DestFileSG)
		DestFileSG.close()
		strdef timefileSG
		timefileSG = "SG_Cell_1_Times.dat"
		TimeFileSG.wopen(timefileSG)
		OutputTimesSG.fwrite(TimeFileSG)
		TimeFileSG.close()
	//}
	ReadFileSG.close()
	objref TempVectorSG
	print "SG Done"
}


objref DestFileZ2, TimeFileZ2, ReadFileZ2, ReadTimeFileZ2, TempVectorZ2 //ADDITION: Time File for recording of spike times out of different SG Cells.
strdef SaveToMeZ2

proc saveZ2(){ local k
	DestFileZ2 = new File()
	ReadFileZ2 = new File()
	TimeFileZ2 = new File() // File for AP Time vector storage.
	TempVectorZ2 = new Vector()
	ReadFileZ2.ropen("SGSCS_Cell_filenames.dat") //MAKE SURE THIS NAME MATCHES YOUR LIST OF OUTPUT FILES.	
	
	//for k=0, 2 { Save multiple neurons. OPTIONAL.  Will need to add bracketed index entry to TempVector (e.g. "TempVector[k]").  
		TempVectorZ2 = OutputVectorsZ2
		ReadFileZ2.scanstr(SaveToMeZ2)
		DestFileZ2.wopen(SaveToMeZ2)
		TempVectorZ2.fwrite(DestFileZ2)
		DestFileZ2.close()
		strdef timefileZ2
		timefileZ2 = "SGSCS_Cell_1_Times.dat"
		TimeFileZ2.wopen(timefileZ2)
		OutputTimesZ2.fwrite(TimeFileZ2)
		TimeFileZ2.close()
	//}
	ReadFileZ2.close()
	objref TempVectorZ2
	print "Z2 Done"
}


objref DestFileT, TimeFileT, ReadFileT, TempVectorT //ADDITION: Time File for recording of spike times out of different T Cells.
strdef SaveToMeT

proc saveT(){ local k
	DestFileT = new File()
	ReadFileT = new File()
	TimeFileT = new File() // File for AP Time vector storage.
	TempVectorT = new Vector()
	ReadFileT.ropen("T_Cell_filenames.dat") //MAKE SURE THIS NAME MATCHES YOUR LIST OF OUTPUT FILES.
	
	//for k=0, 2 {  Save multiple neurons. OPTIONAL.  Will need to add bracketed index entry to TempVector (e.g. "TempVector[k]").  
		TempVectorT = OutputVectorsT
		ReadFileT.scanstr(SaveToMeT)
		DestFileT.wopen(SaveToMeT)
		TempVectorT.fwrite(DestFileT)
		DestFileT.close()
		strdef timefileT
		timefileT = "T_Cell_1_Times.dat"
		TimeFileT.wopen(timefileT)
		OutputTimesT.fwrite(TimeFileT)
		TimeFileT.close()
	//}
	ReadFileT.close()
	objref TempVectorT
	print "T Done"
}


objref DestFileEX, TimeFileEX, ReadFileEX, TempVectorEX //ADDITION: Time File for recording of spike times out of different T Cells.
strdef SaveToMeEX

proc saveEX(){ local k
	DestFileEX = new File()
	ReadFileEX = new File()
	TimeFileEX = new File() // File for AP Time vector storage.
	TempVectorEX = new Vector()
	ReadFileEX.ropen("EX_Cell_filenames.dat") //MAKE SURE THIS NAME MATCHES THAT GENERATED USING FilenameGenerator.m
	//for k=0, 2 {  Only save 3 motor neurons.
		TempVectorEX = OutputVectorsEX
		ReadFileEX.scanstr(SaveToMeEX)
		DestFileEX.wopen(SaveToMeEX)
		TempVectorEX.fwrite(DestFileEX)
		DestFileEX.close()
		strdef timefileEX
		timefileEX = "EX_Cell_1_Times.dat"
		TimeFileEX.wopen(timefileEX)
		OutputTimesEX.fwrite(TimeFileEX)
		TimeFileEX.close()
	//}
	ReadFileEX.close()
	objref TempVectorEX
	print "EX Done"
}