/* DENTATE GYRUS NETWORK TEMPLATE
TODO:
- Consider allowing seed change via argument to init()
- Change the scaling to not be from 0-100 but from 0-1
- Account for different networks in file output names
- Implement PP stimulation of MCs
- Might rename `global_seed`
- Where to put the following?
// for debug purpose only...
objref netcon, nil
objref vec_debug_
objref PP_init_rnd // random number generator to select an PP stimulus site..
objref nslist, rslist, rs_noiselist
objref ns, rs
- Update docstrings on nc_append_rec() and nc_append()
- Clean the following procs/funcs, which were not changed much from Yim et al.:
- onnect_pp_to_gc()
- connect_pp_to_bc()
- connect_pp_to_hipp()
- Address the `KillMC` aspects that are not integrated in Yim, but from Santhakumar
- Functions to be renamed for consistent style:
- nc_append()
- nc_append_rec()
- Variables to be renamed for consistent style:
- PP_box_startid_
- PP_box_nr_
- DGVt_name_
- print_Vtrace
- print_Raster
*/
// Load ranstream
err_ = load_file("ranstream.hoc")
begintemplate DentateGyrus
// Public methods
public init, run
// Simulation related declarations
// Random number generators
objref rdgc2gc, rdgc2bc, rdgc2mc, rdgc2hc
objref rdbc2gc, rdbc2bc, rdbc2mc
objref rdmc2gc1, rdmc2gc2, rdmc2bc, rdmc2mc, rdmc2mc1, rdmc2hc
objref rdhc2gc, rdhc2bc, rdhc2mc
objref rnd_pp2gc, rnd_pp2bc, rnd_pp2hc
objref rdsynb, rdsyna
objref rdsa
// Declare cvode
objref cvode
// Data management related declarations
public VmT, VmMat, idname
objref Spike[726] // Binary representations of spikes
objref Spike_times[726] // Spike times for each cell
objref VmT // Vector storing time index
objref VmMat[726] // One vector per cell storing membrane potential
objref efile // Membrane potential storage file
// Declarations related to debugging
objref vec_stim[100], vec_stim_noise[500]
objref distri_gc_input_
// Network related declarations
public ngcell, nbcell, nmcell, nhcell, npp, nmcell_per_lamella, ntotal
public p_sprouted_, scale_gpas_dg_, scale_sk_dg_, scale_kir_, scale_gabaa_
public scale_PP_strength, scale_PP2BC_strength, scale_PP2HC_strength
public scale_HC2GC_strength, scale_MC2GC_strength, scale_GC2BC_strength
public scale_BC2GC_strength
public cells, nclist
objref cells , Gcell[500], Bcell[6], Mcell[15], Hcell[6]
objref PPSt[100], PPSt_noise[500]
objref rs, rs_noise, rslist, rs_noiselist, pp_start, pp_noise_start
objref nclist, netcon, netcon_d[100], netcon_d_noise[500]
objref vgc2gc, vgc2bc, vgc2mc, vgc2hc
objref vbc2gc, vbc2bc, vbc2mc
objref vmc2gc, vmc2bc, vmc2mc, vmc2hc
objref vhc2gc, vhc2bc, vhc2mc
objref pp2gc, pp2bc, pp2hc
create aacell // artificial cell if windows (boxed) activation [ TODO ] - Study this more
proc init() {
/* Initialization of the DentateGyrus object
Arguments:
$s1 : str : A label for this simulation
$2 : int : Random state (global RNG seed)
$s3 : str : Name of the group being simulated
*/
strdef idname, group
idname = $s1
random_state = $2
group = $s3
// Specify directory to save data
// [ TODO ] - Allow this to be changed optionally
// without having to specify it
// every time. I.e. need to look at
// how one can specify default
// arguments to functions in NEURON
strdef datadir
datadir = "data/dgnetwork/"
// Simulation and data storage parameters
set_simulation_parameters(random_state)
set_data_management_params()
// Set parameters related to network
set_neuron_population_params()
set_perforant_path_input_params()
set_connectivity_params()
// Initialize pseudorandom number generators
set_pseudorandom_number_generators()
// Parameter modifications
oscillation_param_settings(group)
//yim_param_modifications(fig)
}
// ########################################################################
// FUNCTIONS FOR SETTING PARAMETERS, HOUSEKEEPING, ETC
// ########################################################################
proc set_simulation_parameters() {
/* Sets parameters for numerical aspects of the simulation
Arguments:
$1 : int : Global RNG seed
*/
n_patterns = 12 // number of patterns to simulate (should be 12 for the full sim)
dt = 0.1 // simulation time step [ms]
tstop = 2000 // total simulation time [ms]
trial = $1 // trialnumber = seed of simulation [ GENERALIZE ]
trial_noise = 1 // seed for noise to PP
rseed = $1 // pseudorandom number generator seed
debug_ = 2
spontaneous_activity_ = 0 // whether the GC neuron is spontaneously active
test_sa = 0 // 1 if we are testing spontaneous GC activity (cuts all inputs to GCs)
}
proc set_pseudorandom_number_generators(){
/* Create the pseudorandom number generators for the network
TODO:
- Pass argument to ranstream rather than requiring external
*/
// Create RNGs for network connectivity establishment
// GC -> {GC, BC, MC, HC}
rdgc2gc = new Random(rseed)
rdgc2bc = new Random(rseed)
rdgc2mc = new Random(rseed)
rdgc2hc = new Random(rseed)
rdgc2gc.discunif(-50, 50)
rdgc2bc.discunif(-1,1)
rdgc2mc.discunif(0,2)
rdgc2hc.discunif(-2,2)
// BC -> {GC, BC, MC}
rdbc2gc = new Random(rseed)
rdbc2bc = new Random(rseed)
rdbc2mc = new Random(rseed)
rdbc2gc.discunif(-70, 70)
rdbc2bc.discunif(-1, 1)
rdbc2mc.discunif(-3, 3)
// MC -> {GC, BC, MC, HC}
rdmc2gc1 = new Random(rseed)
rdmc2gc2 = new Random(rseed)
rdmc2bc = new Random(rseed)
rdmc2mc = new Random(rseed)
rdmc2mc1 = new Random(rseed)
rdmc2hc = new Random(rseed)
rdmc2gc1.discunif(25, 175)
rdmc2gc2.discunif(-175, -25)
rdmc2bc.discunif(-3,3)
rdmc2mc.discunif(ngcell+nbcell, ngcell+nbcell+nmcell-1)
rdmc2mc1.discunif(-3, 3)
rdmc2hc.discunif(-2, 2)
// HC -> {GC, BC, MC}
rdhc2gc = new Random(rseed)
rdhc2bc = new Random(rseed)
rdhc2mc = new Random(rseed)
rdhc2gc.discunif(-130, 130)
rdhc2bc.discunif(-2, 2)
rdhc2mc.discunif(-2, 2)
// PP -> {GC, BC, HC}
rnd_pp2gc = new Random(rseed)
rnd_pp2bc = new Random(rseed)
rnd_pp2hc = new Random(rseed)
rnd_pp2gc.discunif(0,npp-1)
rnd_pp2bc.discunif(0,npp-1)
rnd_pp2hc.discunif(0,npp-1)
// RNGs for synapses
rdsyna = new Random(rseed)
rdsynb = new Random(rseed)
rdsyna.discunif(0, 1)
rdsynb.discunif(0, 3)
// RNG for spontaneous activity
rdsa = new Random(rseed)
rdsa.discunif(0, ngcell-1)
}
proc set_data_management_params(){
/* Parameters for how data are managed during the simulation.
*/
// Define file names for simulation results output
strdef suffix
suffix = "txt"
// Flags about what to print
print_Vtrace = 0 // print voltage trace to file
print_Raster = 1 // print spike raster to file
print_template = 1 // write out spike sequence into template file
print_GC_sprouting_input_ = 0 // calculates and print out GC sprouting input number distribution
print_stim_ = 1 // print stimulus to GCs for debug purpose..
print_stim_noise = 1 // print noise to GCs for debug purpose..
}
proc set_neuron_population_params(){
/* Initialize the number of cells of each type in the networks
TODO:
- Generalize to allow these to vary
*/
ngcell = 500 // number of GCs
nbcell = 6 // number of BCs (Samathakumar 2005)
nmcell = 15 // number of MCs
nhcell = 6 // number of HCs (Samathakumar 2005)
npp = 100 // ECII Neurons (Myers and Scharfman, 2009)
nmcell_per_lamella = 3 // Number of mossy cells per lamella
// Add 2*npp because of PP input and PP noise...
ntotal = ngcell + nbcell + nmcell + nhcell + npp + ngcell
// Neuronal params
v_init = -77.71
}
proc set_perforant_path_input_params() {
/* Set parameters for the perforant path stimulation
*/
PP_nstimulus_ = 1 // one input per GC
PP_input2MC_ = 0 // stimulate MC with PP
PP_rate_ = 10. // rate of PP Poisson input
PP_rate_noise_ = 0. // rate of PP Poisson noise
PP_box_nr_ = 6 // number of active PPs
PP_box_stop_ = 2000. // time of box [ms]
PP_box_start_ = 1. // shift of box [ms]
PP_box_nspk_ = 25 // number of spike per active PP
PP_freq_ = 3 // oscillation frequency
if (PP_input2MC_ != 0) {
print "PP stimulation not yet implemented if stimulation of also MCs...."
quit()
}
}
proc set_connectivity_params() {
/* Parameters related to network connectivity
Arguments:
*/
scale_factor = 1.75 // coefficient for scale PP and scale PP2BC and PP2GC. From rate model analysis.
p_sprouted_ = 0 // proportion of other GCs to which a GC will connect
scale_gpas_dg_ = 1 // scaling [%/100] of gpas of th GC model
scale_sk_dg_ = 1 // scaling of Ca-dependent K (SK) of th GC model
scale_kir_ = 1 // scaling of KIR conductance in GC model
scale_gabaa_ = 1 // scaling of GABAA in GC model
scale_PP_strength = scale_factor*(1-test_sa) // scaling of synaptic weight of PP
scale_PP2BC_strength = 1*(1-test_sa) // scaling of synaptic weight PP->BC (100% for original value which was 50% of PP->GC strength..)
scale_PP2HC_strength = 0 // scaling of synaptic weight PP->HIPP (set to 0% because there are no PP->HIPP connections...)
scale_HC2GC_strength = 1*(1-test_sa) // scaling of synaptic weight HC->GC (beta_HIPP in Myers and Scharfman, 2009)
scale_MC2GC_strength = 1*(1-test_sa) // scaling of synaptic weight MC->GC (beta_HIPP in Myers and Scharfman, 2009)
scale_GC2BC_strength = 3*(1-test_sa) // scaling of synaptic weight GC->BC
scale_BC2GC_strength = 3*(1-test_sa) // scaling of synaptic weight BC->GC
spontaneous_activity_strength = 0 // scaling of the strength of synaptic weight for the spontaneous AP generators
}
proc yim_param_modifications() {
/* Reset parameters to generate the figures in Yim et al (2015)
Arguments:
$1 : int : Which figure you are trying to generate data for
*/
fig = $1
if (fig == 1) {
scale_PP_strength = 0.10
scale_kir_ = 1
scale_gabaa_ = 1
p_sprouted_ = 0
}
if (fig == 2) {
scale_PP_strength = 0.1
scale_kir_ = 1
scale_gabaa_ = 1
p_sprouted_ = 0.3
}
if (fig == 3) {
scale_PP_strength = 0.1
scale_kir_ = 4.
scale_gabaa_ = 4.
p_sprouted_ = 0.3
}
if (fig == 4) {
scale_PP_strength = 0.16
scale_kir_ = 1.
scale_gabaa_ = 1.
p_sprouted_ = 0.1
}
if (fig == 5) {
scale_PP_strength = 0.16
scale_kir_ = 4.
scale_gabaa_ = 4.
p_sprouted_ = 0.1
}
}
proc oscillation_param_settings() {
/* Sets different parameters for oscillation related
experiments
Arguments:
$s1 : str : Which frequency we are simulating
*/
// Theta oscillations
if (strcmp($s1, "theta") == 0) {
PP_freq_ = 3 // Zheng et al.
PP_box_nspk_ = 70
PP_min_invl_ = 5 // 10 optimal for theta
PP_scale_max_invl_ = 100
scale_na_conductances = 2.5 // gives us U-shaped curve
scale_kdr_conductances = 2.5 // gives us U-shaped curve
scale_ka_conductances = 1
gbar_ht_ = 0.0004
gbar_lt_ = 0
scale_size_ = 1
scale_gpas_dg_ = 1
scale_sk_dg_ = 1
scale_gabaa_ = 1
scale_kir_ = 0
spontaneous_activity_rate = 0.25 // In Hz
}
// Fast theta oscillations
if (strcmp($s1, "ftheta") == 0) {
PP_freq_ = 8 // Zheng et al.
PP_box_nspk_ = 100
PP_min_invl_ = 5 // 10 optimal for theta
PP_scale_max_invl_ = 100
scale_na_conductances = 2.5 // gives us U-shaped curve
scale_kdr_conductances = 2.5 // gives us U-shaped curve
scale_ka_conductances = 1
gbar_ht_ = 0.0004
gbar_lt_ = 0
scale_size_ = 1
scale_gpas_dg_ = 1
scale_sk_dg_ = 1
scale_gabaa_ = 1
scale_kir_ = 0
spontaneous_activity_rate = 0.25 // In Hz
}
// Alpha oscillations
if (strcmp($s1, "alpha") == 0) {
PP_freq_ = 12 // Zheng et al.
PP_box_nspk_ = 120
PP_min_invl_ = 5
PP_scale_max_invl_ = 100
scale_na_conductances = 2.5 // gives us U-shaped curve
scale_kdr_conductances = 2.5 // gives us U-shaped curve
scale_ka_conductances = 1
gbar_ht_ = 0.0004
gbar_lt_ = 0
scale_size_ = 1
scale_gpas_dg_ = 1
scale_sk_dg_ = 1
scale_gabaa_ = 1
scale_kir_ = 0
spontaneous_activity_rate = 0.25 // In Hz
}
// Beta oscillations
if (strcmp($s1, "beta") == 0) {
PP_freq_ = 20 // Zheng et al.
PP_box_nspk_ = 200
PP_min_invl_ = 5
PP_scale_max_invl_ = 100
scale_na_conductances = 2.5 // gives us U-shaped curve
scale_kdr_conductances = 2.5 // gives us U-shaped curve
scale_ka_conductances = 1
gbar_ht_ = 0.0004
gbar_lt_ = 0
scale_size_ = 1
scale_gpas_dg_ = 1
scale_sk_dg_ = 1
scale_gabaa_ = 1
scale_kir_ = 0
spontaneous_activity_rate = 0.25 // In Hz
}
// Gamma oscillations
if (strcmp($s1, "gamma") == 0) {
PP_freq_ = 35 // Zheng et al.
PP_box_nspk_ = 300
PP_min_invl_ = 3 // 1 optimal for gamma
PP_scale_max_invl_ = 100
scale_na_conductances = 2.5 // gives us U-shaped curve
scale_kdr_conductances = 2.5 // gives us U-shaped curve
scale_ka_conductances = 1
gbar_ht_ = 0.0004
gbar_lt_ = 0
scale_size_ = 1
scale_gpas_dg_ = 1
scale_sk_dg_ = 1
scale_gabaa_ = 1
scale_kir_ = 0
spontaneous_activity_rate = 0.25 // In Hz
}
}
// ########################################################################
// FUNCTIONS FOR SETTING UP CELLS/INPUTS
// ########################################################################
proc make_cells() {local i
// Need to re-declare cell arrays with proper sizes
objref Gcell[ngcell], Bcell[nbcell], Mcell[nmcell], Hcell[nhcell]
objref PPSt[npp], PPSt_noise[ngcell]
cells = new List()
// Create populations
for i=0, ngcell-1 {Gcell[i] = new GranuleCell(i,scale_na_conductances, scale_kdr_conductances, scale_ka_conductances, gbar_ht_, gbar_lt_, scale_size_, scale_gpas_dg_, scale_sk_dg_, scale_gabaa_, scale_kir_)}
for i=0, nbcell-1 {Bcell[i] = new BasketCell(i)}
for i=0, nmcell-1 {Mcell[i] = new MossyCell(i)}
for i=0, nhcell-1 {Hcell[i] = new HIPPCell(i)}
for i=0, npp-1 {PPSt[i] = new PPstim(i, PP_rate_, tstop, PP_box_start_, PP_box_stop_, PP_freq_, PP_min_invl_, PP_scale_max_invl_)}
// [TESTING] Adding individual stims to each cell in order to simulate spontaneous activity
//for i=0, npp-1 {PPSt_noise[i] = new PPstim(i, PP_rate_, tstop, PP_box_start_, PP_box_stop_)}
for i=0, ngcell-1 { PPSt_noise[i] = new PPstim(i, PP_rate_, tstop, 1., tstop-1, 0, 10, 10) }
// Append to cells list
for i = 0, ngcell-1 {cells.append(Gcell[i])} // cells 0-499 GCs
for i = 0, nbcell-1 {cells.append(Bcell[i])} // cells 500-505 BC
for i = 0, nmcell-1 {cells.append(Mcell[i])} // cells 506-520 MC
for i = 0, nhcell-1 {cells.append(Hcell[i])} // cells 521-526 HC
for i = 0, npp-1 {cells.append(PPSt[i])} // 527 - xxx PP artificial cell
// [TESTING] Append GC spontaneous activity generators to list
//for i = 0, npp-1 {cells.append(PPSt_noise[i])}
if (spontaneous_activity_ == 1) {
for i = 0, ngcell-1 { cells.append(PPSt_noise[i]) }
}
}
proc make_input_stimulation() { local i
/* Creates the perforant path input stimulation objects
*/
// Print update to console
print "\tDefining Random Generator for each PPStim object."
objref pp_start, pp_noise_start
objref netcon_d[npp]
objref netcon_d_noise[ngcell]
rslist = new List()
rs_noiselist = new List()
// Input window (box)
random_stream_offset_ = PP_rate_ * 1000
for i=0, npp-1 {
rs = new RandomStream(i, random_state)
rslist.append(rs)
PPSt[i].pp.noiseFromRandom(rs.r)
rs.r.uniform(0,1)
rs.start(random_state)
}
// artificial pulse to PPSt[i].pp in order to become active...
aacell pp_start = new NetStim125(.5)
pp_start.interval = 1e-5
pp_start.number = 1
pp_start.start = 0
pp_start.forcestop = 1.
pp_start.noise = 1
// assign Random Generator to init..
// [TODO] - Study more what is happening here
rs = new RandomStream(npp, random_state)
return_val_ = rslist.append(rs) // save returned value in return_val_ to suppress screen output
return_val_ = pp_start.noiseFromRandom(rs.r)
return_val_ = rs.r.negexp(1)
return_val_ = rs.start(random_state)
// each PPSt[i].pp needs to receive a netevent to become active...
// [ TODO ] - allow for these parameters to be changed more generally
for i=0, npp-1 {
netcon_d[i] = new NetCon(pp_start,PPSt[i].pp)
netcon_d[i].weight = 10.
netcon_d[i].delay = 0.001
netcon_d[i].threshold = 10.
}
if (spontaneous_activity_ == 1) {
// Input window (box)
for i=0, ngcell-1 {
rs_noise = new RandomStream(i, random_state)
rs_noiselist.append(rs_noise)
PPSt_noise[i].pp.noiseFromRandom(rs_noise.r)
rs_noise.r.uniform(0,1)
rs_noise.start(random_state)
}
// artificial pulse to PPSt_noise[i].pp in order to become active...
aacell pp_noise_start = new NetStim125(.5)
pp_noise_start.interval = 1e-5
pp_noise_start.number = 1
pp_noise_start.start = 0
pp_noise_start.forcestop = 1.
pp_noise_start.noise = 1
// assign random number generator
rs_noise = new RandomStream(ngcell, random_state)
return_val_ = rslist.append(rs_noise) // save returned value in return_val_ to suppress screen output
return_val_ = pp_noise_start.noiseFromRandom(rs_noise.r)
return_val_ = rs_noise.r.negexp(1)
return_val_ = rs_noise.start(random_state)
// each PPSt_noise[i].pp needs to receive a netevent to become active...
for i=0, ngcell-1 {
netcon_d_noise[i] = new NetCon(pp_noise_start,PPSt_noise[i].pp)
netcon_d_noise[i].weight = 10.
netcon_d_noise[i].delay = 0.001
netcon_d_noise[i].threshold = 10.
}
}
}
// ########################################################################
// FUNCTIONS FOR MAKING NETWORK CONNECTIVITY
// ########################################################################
proc make_connections(){
/* Instantiates connectivity in the network
*/
// Instantiate the connectivity list
nclist = new List()
// Instantiate vectors for connections
// GC -> {GC, BC, MC, HC}
vgc2gc = new Vector(ngcell, 0)
vgc2bc = new Vector(nbcell, 0)
vgc2mc = new Vector(nmcell, 0)
vgc2hc = new Vector(nhcell, 0)
// BC -> {GC, BC, MC}
vbc2gc = new Vector(ngcell, 0)
vbc2bc = new Vector(nbcell, 0)
vbc2mc = new Vector(nmcell, 0)
// MC -> {GC, BC, MC, HC}
vmc2gc = new Vector(ngcell, 0)
vmc2bc = new Vector(nbcell, 0)
vmc2mc = new Vector(nmcell, 0)
vmc2hc = new Vector(nhcell, 0)
// HC -> {GC, BC, MC}
vhc2gc = new Vector(ngcell, 0)
vhc2bc = new Vector(nbcell, 0)
vhc2mc = new Vector(nmcell, 0)
// Connect areas
connect_pp_to_gc()
connect_pp_to_bc()
connect_pp_to_hipp()
connect_granule_cells()
connect_basket_cells()
connect_mossy_cells()
connect_hipp_cells()
if (spontaneous_activity_ == 1) {
add_spontaneous_gc_activity()
}
}
func nc_append() {
// neuron connect $1 with $2.pre_list.object($3), weight $4, delay $5, threshold $6
// connects:
// cells.object($1) with
// $o1 = cells.object($2).pre_list.object($3) and
// returns:
// netcon = $o2
if ($3 >= 0 ) {
// connect_pre is function in the respective cell definition
cells.object($1).connect_pre(cells.object($2).pre_list.object($3),netcon)
netcon.weight = $4
netcon.delay = $5
netcon.threshold = $6
}
nclist.append(netcon)
return nclist.count-1
}
func nc_append_rec() {
/* neuron connect $1 with $2.pre_list.object($3), weight $4, delay $5, threshold $6
// connects:
// cells.object($1) with
// $o1 = cells.object($2).pre_list.object($3) and
// returns:
// netcon = $o2
// record events to $o7
*/
if ($3 >= 0 ) {
// connect_pre is function in the respective cell definition
cells.object($1).connect_pre(cells.object($2).pre_list.object($3),netcon)
netcon.weight = $4
netcon.delay = $5
netcon.threshold = $6
netcon.record($o7)
}
nclist.append(netcon)
return nclist.count-1
}
func is_connected() {local i, c localobj net_c
/* Checks for preexisting connections between randomly selected cells
to avoid multiple contacts between same 2 cells
*/
c=0
for i=0, nclist.count-1 {
net_c = nclist.object(i)
if (($o1 == net_c.postcell()) && ($o2 == net_c.precell())) {c=1}
}
return c
}
proc connect_pp_to_gc() { local i, j localobj pprec
/* Function that connects perforant path inputs to granule cells.
*/
// Print message to indicate procedure
print "Connecting PP to post-synaptic targets."
print "\tConnecting PP -> GC."
// Create vector that marks PP neurons that already project to a GC
pprec = new Vector(npp, 0)
// Objects to record inputs/noise for debugging [TODO] verify necessity...
// Must be re-declared with proper value of npp
objref vec_stim[npp]
for i = 0, npp-1 { vec_stim[i] = new Vector() }
// Make connections for each individual GC
for i=0, ngcell-1 {
// Re-initialize connectivity vector to prevent multiple stimulation
// connections from the same PP cell
pp2gc = new Vector(npp, 0)
// Each GC receives input from random set of 20% of PP neurons
// [ TODO ] - ALLOW TO HAVE 20% CHANGED TO DIFFERENT VALUES
for nr_conn = 0, int(npp/5.) - 1 {
// Pick a random PP cell to connect to the GC
// Repeat until j is a PP cell that is not yet connected
// Then mark j as connected so it is not chosen again
j = rnd_pp2gc.repick()
while (pp2gc.x[j] == 1) {
j = rnd_pp2gc.repick()
}
pp2gc.x[j] = 1
if ((print_stim_==1) && (i<npp)) {
if (pprec.x[j] == 0) {
// record input sequence to first neuron for debug purpose only....
nc_append_rec(j+ngcell+nbcell+nmcell+nhcell, i, 0, 2e-2*scale_PP_strength, 3, 10, vec_stim[j])
pprec.x[j] = 1
} else {
// record input sequence to first neuron for debug purpose only....
nc_append(j+ngcell+nbcell+nmcell+nhcell, i, 0, 2e-2*scale_PP_strength, 3, 10)
}
} else {
// connect PP to GC[j],syn[0],wt,del,threshold <AH> NOTE both synapses are equal, ie. weight delay (except position) not important
nc_append(j+ngcell+nbcell+nmcell+nhcell, i, 0, 2e-2*scale_PP_strength, 3, 10)
}
}
// Print debugging message
if (debug_ == 3) {
print "\nGC: ",i
for ii=0,pp2gc.size()-1 { printf ("%d, ",pp2gc.x[ii]) }
}
}
}
proc add_spontaneous_gc_activity() { local i,j
/* Function that adds spontaneous hyperactivity to the granule cells
*/
// Print message to indicate procedure
print "Adding spontaneous GC activity."
// Objects to record inputs/noise for debugging [TODO] verify necessity...
// Must be re-declared with proper value of npp
objref vec_stim_noise[ngcell]
for i = 0, ngcell-1 { vec_stim_noise[i] = new Vector() }
// Add spontaneous activity generator for each individual GC
for i=0, ngcell-1 {
if ((print_stim_ == 1)) {
nc_append_rec(i+ngcell+nbcell+nmcell+nhcell+npp, i, 0, 2e-2*spontaneous_activity_strength, 0, 0, vec_stim_noise[i])
} else {
nc_append(i+ngcell+nbcell+nmcell+nhcell+npp, i, 0, 2e-2*spontaneous_activity_strength, 0, 0)
}
}
}
proc connect_pp_to_bc() { local i,j
/* Connect perforant path inputs to basket cells
*/
// Print message to indicate procedure
print "\tConnecting PP -> BC."
// Only proceed if there are PP to BC connections
if (scale_PP2BC_strength != 0) {
for i=ngcell, ngcell+nbcell-1 {
pp2bc = new Vector(npp, 0)
// Each BC receives input from random set of 20% of PP neurons
// [ TODO ] - ALLOW TO HAVE 20% CHANGED TO DIFFERENT VALUES
for nr_conn = 0, int(npp/5.)-1 {
j = rnd_pp2bc.repick()
while (pp2bc.x[j] == 1) { j = rnd_pp2bc.repick() }
pp2bc.x[j] = 1
// connect PP to BC[j],syn[0],wt,del,threshold <AH> NOTE both synapses are equal
nc_append(j+ngcell+nbcell+nmcell+nhcell, i, 0, 1e-2*scale_PP_strength*scale_PP2BC_strength, 3, 10)
}
// Print debugging message
if (debug_ == 3) {
print "\nBC: ",i
for ii=0,pp2bc.size()-1{ printf ("%d, ",pp2gc.x[ii])}
}
}
}
}
proc connect_pp_to_hipp() { local i,j
/* Connecting perorant path to HIPP cells
*/
if (scale_PP2HC_strength != 0) {
print "\tConnecting PP -> HIPP."
for i=0, nhcell-1 {
pp2hc = new Vector(npp, 0)
// Each HIPP receives input from random set of 20% of PP neurons
// [ TODO ] - ALLOW TO HAVE 20% CHANGED TO DIFFERENT VALUES
for nr_conn = 0, int(npp/5.)-1 {
j = rnd_pp2hc.repick()
while (pp2hc.x[j] == 1) { j = rnd_pp2hc.repick() }
pp2hc.x[j] = 1
// used here synaptic strength of GC->HIPP / connections PP-> HIPP not in Santhakumar 2005
nc_append(j+ngcell+nbcell+nmcell+nhcell, i+ngcell+nbcell+nmcell, 0, 0.5e-3*scale_PP_strength*scale_PP2HC_strength, 3, 10)
}
if (debug_ == 3) {
print "\nHIPP: ",i
for ii=0,pp2bc.size()-1{ printf ("%d, ",pp2hc.x[ii])}
}
}
if (debug_ == 3) { print "PP -> HIPP -DONE-" }
}
}
proc connect_granule_cells() { local i,j localobj c1, c2
/* Connect granule cells to post-synaptic targets
*/
// Print to indicate step
print "Connecting GCs to post-synaptic targets."
if (debug_ == 0) { print "n_sprout_ = ",n_sprout_ }
if (p_sprouted_>0 && debug_ == 2) {print "\tNOTE: Implementing mossy fibre sprouting."}
for i=0, ngcell-1 {
// CONNECT GC'S TO BC'S
for j=0, 0 {
// Which lamella is GC[i] in wrt basket cells
// Note: one basket cell per lamella
// [ TODO ] - Allow for more than one BC per lamella
if (i < ngcell/nbcell) { a=0 }
if ((i > ngcell/nbcell-1) && (i < ngcell*2/nbcell)) { a=1 }
if ((i > ngcell*2/nbcell-1) && (i < ngcell*3/nbcell)) { a=2 }
if ((i > ngcell*3/nbcell-1) && (i < ngcell*4/nbcell)) { a=3 }
if ((i > ngcell*4/nbcell-1) && (i < ngcell*5/nbcell)) { a=4 }
if ((i > ngcell*5/nbcell-1) && (i < ngcell)) { a=5}
// Randomly pick location of post synaptic Bcell from distribution [-1:1]
Gauz3 = rdgc2bc.repick()
// Determine appropriate post syn BC
if (a+Gauz3 > nbcell-1) { npost = a+Gauz3-nbcell }
if (a+Gauz3 < 0) { npost = a+Gauz3+nbcell }
if ((a+Gauz3 > -1) && (a+Gauz3 < nbcell)) {npost = a+Gauz3}
// randomly pick the dendrite to connect to from [0:3] ( i.e. // randomize among 4 dendrites )
dbr = rdsynb.repick()
if (debug_ == 1 ) {
print "GC \t",i,"\t to BC\t\t",npost, a
}
// check to make sure that post syn BC does not get more than 90 GC synapses
// [ TODO ] - ALLOW FOR DIFFERENT NUMBER OF MAX SYNAPSES
if (vgc2bc.x[npost] < 90) {
// connect GC[i] to BC[j],syn[2]+dendritic_var,wt,del,threshold
nc_append(i, ngcell+npost, dbr+2, 4.7e-3*scale_GC2BC_strength, .8, 10)
if (debug_ == 0 ) {
print "GC \t",i,"\t to BC\t\t",npost,"\t",dbr+2
}
// increment the no of synapses to the post cell
vgc2bc.x[npost] +=1
} else {
j -= 1
if (debug_ == 1 ) {print "nogc2bc"}
} // for connection that is not made reconnect axon to another cell
}
// CONNECT GC'S TO MC'S
for j=0, 0 {
// Based on the lamellar distribution of the GCs to MCs - 500 GCs were divided into 5 groups, 3 MCs were distributed in each lamella
// print "Based on the lamellar distribution of the GCs to MCs..."
// [ TODO ] - Allow for different number of MC's per lamella
if (i < ngcell/5) { a=0}
if ((i > ngcell/5-1) && (i < ngcell*2/5)) { a=1}
if ((i > ngcell*2/5-1) && (i < ngcell*3/5)) { a=2}
if ((i > ngcell*3/5-1) && (i < ngcell*4/5)) { a=3}
if ((i > ngcell*4/5-1) && (i < ngcell)) { a=4}
b=a*3 // from [0:12]
npost = rdgc2mc.repick() // from [0:2]
dbr = rdsynb.repick() // from [0:2]
// [ TODO ] - ALLOW FOR DIFFERENT NUMBER OF MAX SYNAPSES
if (vgc2mc.x[npost+b] < 38){
nc_append(i, ngcell+nbcell+npost+b, dbr+4, 0.2e-3, 1.5, 10)
if (debug_ == 1 ) {
print "GC \t",i,"\t to MC\t\t", npost+b, "\t", dbr+4
}
vgc2mc.x[npost+b] +=1
} else {
j -= 1
if (debug_ == 1 ) {
print "nogc2mc"
}
}
}
// CONNECT GC's to HIPP
for j=0, 2 {
// <AH> comment added:
// Based on the lamellar distribution of the GCs to HIPPs - 500 GCs were divided into 6 groups
if (i < ngcell/6) { a=0}
if ((i > ngcell/6-1) && (i < ngcell*2/6)) { a=1}
if ((i > ngcell*2/6-1) && (i < ngcell*3/6)) { a=2}
if ((i > ngcell*3/6-1) && (i < ngcell*4/6)) { a=3}
if ((i > ngcell*4/6-1) && (i < ngcell*5/6)) { a=4}
if ((i > ngcell*5/6-1) && (i < ngcell)) { a=5}
Gauz3 = rdgc2hc.repick()
if (a+Gauz3 > nhcell-1) { npost = a+Gauz3-nhcell }
if (a+Gauz3 < 0) { npost = a+Gauz3+nhcell }
if ((a+Gauz3 > -1) && (a+Gauz3 < nhcell)) { npost = a+Gauz3 }
dbr = rdsynb.repick()
c1 = cells.object(ngcell+nbcell+nmcell+npost)
c2 = cells.object(i)
if ((is_connected(c1, c2) == 0) && (vgc2hc.x[npost] < 275)) {
nc_append(i, ngcell+nbcell+nmcell+npost, dbr, 0.5e-3, 1.5, 10)
if (debug_ == 1 ) {
print "GC \t",i,"\t to HC\t\t",npost, "\t", dbr
}
vgc2hc.x[npost] +=1
} else {
j -= 1
}
}
// GCs -> GCs
// NOTE: THIS IS FOR SPROUTED SYNAPSES
// NOTE: 100% Sprouting = 100 new synapses! (compare p. 443 in Santhakumar paper)
n_sprout_ = p_sprouted_ - 1
// 9 in original file // each GC is recurrent connected to 10 GC
// (i.e. 10/500 => p = 0.02) but sprouting is diff different -> see above
for j=0, n_sprout_ {
Gauz3 = rdgc2gc.repick()
if (i+Gauz3 > 499) { npost = i+Gauz3-500 }
if (i+Gauz3 < 0) { npost = i+Gauz3+500 }
if ((i+Gauz3 > -1) && (i+Gauz3 < 500)) { npost = i+Gauz3 }
dbr = rdsyna.repick() // [0:1]
c1 = cells.object(npost)
c2 = cells.object(i)
if ((is_connected(c1, c2) == 0) && (vgc2gc.x[npost] < (n_sprout_*1.5+2) )) { // if is connected AND not more than 14 incoming connections...
// (original file < 15) (assume 1.5 times average connections for upper limit...)
nc_append(i, npost, dbr+7, 2e-3, .8, 10) // Gcell[3] to Bcell[1]
// print npost, dbr+8
if (debug_ == 0 ) {print "GC \t",i,"\t to GC\t\t",npost, "\t", dbr+7}
vgc2gc.x[npost] +=1
} else {
j -= 1
if (debug_ == 0) {print "gc2gc"}
}
}
}
if (print_GC_sprouting_input_ == 1) {
distri_gc_input_ = new Vector()
max_gc_input_ = 0
if (debug_ ==2) { print "Calculate GC-GC Input Distribution"}
for zz = 0, int(n_sprout_*1.5+2) {distri_gc_input_.append(0)}
for npost=0,ngcell-1 {
distri_gc_input_.x[vgc2gc.x[npost]]+=1
if (vgc2gc.x[npost]>max_gc_input_) { max_gc_input_ = vgc2gc.x[npost]} // find max input number
}
for zz = 0, int(n_sprout_*1.5+2) {print zz,"\t",distri_gc_input_.x[zz]}
print "maximum input number is:\t",max_gc_input_
}
}
proc connect_basket_cells() { local i,j localobj c1, c2
print "Connecting BCs to post-synaptic targets. "
for i=0, nbcell-1 {
// BC -> GC
// [ TODO ] - Allow for different lamella sizes
for j=0, 99 {
Gauz3 = rdbc2gc.repick() // [-70:70]
if (debug_ == 1 ) {print Gauz3}
if (i*83+41+Gauz3 > 499) {npost = i*83+41+Gauz3-500 }
if (i*83+41+Gauz3 < 0) {npost = i*83+41+Gauz3+500}
if ((i*83+41+Gauz3 > -1) && (i*83+41+Gauz3 < 500)) {npost = i*83+41+Gauz3}
if (debug_ == 1 ) {print i, npost}
if (nbcell != 6) {max_conn = 4} else {max_conn = 2} // if not original setup use more spread...
c1 = cells.object(npost)
c2 = cells.object(i+ngcell)
if ((is_connected(c1, c2) == 0) && (vbc2gc.x[npost] < max_conn)) { // change < 2 to < 4
nc_append(i+ngcell, npost, 6, 1.6e-3*scale_BC2GC_strength, .85, -10)
vbc2gc.x[npost] +=1
if (debug_ == 1 ) { print i, npost, 6 }
} else {
j -= 1
if (debug_ == 1 ) {print "BC2GC"}
}
}
// BC -> BC
for j=0, 1 {
Gauz3 = rdbc2bc.repick() // [-1,0,1] (postsyn spread around single id...)
//print Gauz3
if (i+Gauz3 > nbcell-1) {npost = i+Gauz3-nbcell } // periodic boundary conditions
if (i+Gauz3 < 0) {npost = i+Gauz3+nbcell}
if ((i+Gauz3 >-1) && (i+Gauz3 < nbcell)) {npost = i+Gauz3}
dbr = rdsyna.repick() // [0:1]
if (nbcell != 6) {max_conn = 4} else {max_conn = 3}
c1 = cells.object(npost+ngcell)
c2 = cells.object(i+ngcell)
if ((is_connected(c1, c2) == 0) && (vbc2bc.x[npost] < max_conn)) { // change < 3 to < 4
nc_append(i+ngcell, npost+ngcell, dbr+8, 7.6e-3, .8, -10)
if (debug_ == 1 ) {print npost, dbr+8}
vbc2bc.x[npost] +=1
} else {
j -= 1
if (debug_ == 1 ) {print "bc2bc"}
}
}
// BC -> MC
for j=0, 2 {
Gauz3 = rdbc2mc.repick() // [-3:3]
if (i*2+2+Gauz3 > 14) {npost = i*2+2+Gauz3-15 }
if (i*2+2+Gauz3 < 0) {npost = i*2+2+Gauz3+15}
if ((i*2+2+Gauz3 >-1) && (i*2+2+Gauz3 < 15)) {npost = i*2+2+Gauz3}
// if ((is_connected(MossyCell[npost], BasketCell[i]) == 0) && (vbc2mc.x[npost] < 3) && (killMC.contains(ngcell+nbcell+npost) == 0)) { // use if killing MC
if (nbcell != 6) {max_conn = 4} else {max_conn = 3}
c1 = cells.object(npost+ngcell+nbcell)
c2 = cells.object(i+ngcell)
if ((is_connected(c1, c2) == 0) && (vbc2mc.x[npost] < max_conn)) { // change < 3 to < 4
nc_append(i+ngcell, npost+ngcell+nbcell, 12, 1.5e-3, 1.5, -10) // Gcell[3] to Bcell[1]
if (debug_ == 1 ) {print npost, 12}
vbc2mc.x[npost] +=1
} else {
j -= 1
if (debug_ == 1 ) {print "bc2mc"}
}
// if (killMC.contains(ngcell+nbcell+npost) == 1) {j +=1 if (debug_ == 1 ) {print "dead MC"}} // use if killing MC
}
}
}
proc connect_mossy_cells() {local i,j localobj c1, c2
/* Connect mossy cells to post-synaptic targets
*/
// Print to indicate step
print "Connecting MCs to post-synaptic targets."
for i=0, nmcell-1 {
//if (killMC.contains(ngcell+nbcell+i) == 0) // use if killing MC
if (i < 3) { y=0 }
if ((i > 2) && (i < 6)) { y=1 }
if ((i > 5) && (i < 9)) { y=2 }
if ((i > 8) && (i < 12)) { y=3 }
if ((i > 11) && (i < 15)) { y=4 }
// MC -> GC1
for j=0, 99 {
Gauz1 = rdmc2gc1.repick() // [25:175]
if (i*33+17+Gauz1 > 499) {
npost1 = i*33+17+Gauz1-500
} else {
npost1 =i*33+17+Gauz1
}
dbr = rdsyna.repick() // [0:1]
c1 = cells.object(npost1)
c2 = cells.object(i+ngcell+nbcell)
if ((is_connected(c1, c2) == 0) && (vmc2gc.x[npost1] < 7)) {
nc_append(i+ngcell+nbcell, npost1, dbr+2, 0.3e-3*scale_MC2GC_strength, 3, 10)
vmc2gc.x[npost1] +=1
} else {
j -= 1
}
}
// MC -> GC2
for j=0, 99 {
Gauz2 = rdmc2gc2.repick() // [-175:25]
if (i*33+17+Gauz2 < 0) {
npost2 =i*33+17+Gauz2+500
} else {npost2 =i*33+17+Gauz2}
dbr = rdsyna.repick() // [0:1]
c1 = cells.object(npost2)
c2 = cells.object(i+ngcell+nbcell)
if ((is_connected(c1, c2) == 0) && (vmc2gc.x[npost2] < 7)) {
nc_append(i+ngcell+nbcell, npost2, dbr+2, 0.3e-3*scale_MC2GC_strength, 3, 10) // Gcell[3] to Bcell[1]
vmc2gc.x[npost2] +=1
} else {
j -= 1
}
}
// MC -> BC
for j=0, 0 {
Gauz3 = rdmc2bc.repick() // Gauz3 = [-3:3]
if (y+Gauz3 > nbcell-1) {npost = y+Gauz3-nbcell} // y = [0:4] => y+Gaus3 = [-3:7]
if (y+Gauz3 < 0) {npost = y+Gauz3+nbcell}
if ((y+Gauz3 > -1) && (y+Gauz3 < nbcell)) {npost = y+Gauz3}
dbr = rdsyna.repick()
if ((vmc2bc.x[npost] < 4) && (Gauz3 !=0)) {
nc_append(i+ngcell+nbcell, ngcell+npost, dbr+6, 0.3e-3, 3, 10) // Gcell[3] to Bcell[1]
vmc2bc.x[npost] += 1
} else {
j -= 1
}
}
// MC -> MC
for j=0, 2 {
Gauz3 = rdmc2mc1.repick() //[-3:3]
if (i+Gauz3 > 14) {npost = i+Gauz3-15 }
if (i+Gauz3 < 0) {npost = i+Gauz3+15}
if ((i+Gauz3 >-1) && (i+Gauz3 < 15)) {npost = i+Gauz3}
dbr = rdsynb.repick()
// if ((is_connected(MossyCell[npost], MossyCell[i]) == 0) && (vmc2mc.x[npost] < 4) && (Gauz3 != 0) && (killMC.contains(ngcell+nbcell+npost) == 0)) // use if killing MC
c1 = cells.object(npost+ngcell+nbcell)
c2 = cells.object(i+ngcell+nbcell)
if ((is_connected(c1, c2) == 0) && (vmc2mc.x[npost] < 4) && (Gauz3 != 0)) {
nc_append(i+ngcell+nbcell, npost+ngcell+nbcell, dbr+8, 0.5e-3, 2, 10)
// print npost, dbr+8
vmc2mc.x[npost] +=1
} else {
j -= 1
}
// if (killMC.contains(ngcell+nbcell+npost) == 1){ j += 1 print "dead MC"} // use if killing MC
}
// MC -> HC
for j=0, 1 {
Gauz3 = rdmc2hc.repick() // [-2:2]
if (y+Gauz3 > nhcell-1) {npost = y+Gauz3-nhcell} // changed code here to allow > 6 HCells
if (y+Gauz3 < 0) {npost = y+Gauz3+nhcell}
if ((y+Gauz3 > -1) && (y+Gauz3 < nhcell)) {npost = y+Gauz3}
dbr = rdsynb.repick()
c1 = cells.object(ngcell+nbcell+nmcell+npost)
c2 = cells.object(i+ngcell+nbcell)
if ((is_connected(c1, c2) == 0) && (vmc2hc.x[npost] < 7) && (Gauz3 != 0)) {
nc_append(i+ngcell+nbcell, ngcell+nbcell+nmcell+npost, dbr+4, 0.2e-3, 3, 10) // Gcell[3] to Bcell[1]
// print npost, dbr+4
vmc2hc.x[npost] +=1
} else {
j -= 1
}
}
}
}
proc connect_hipp_cells() {local i,j localobj c1, c2
/* Connect HIPP cells to post-synaptic targets
*/
// Print to indicate step
print "Connecting HIPP cells to post-synaptic targets."
for i=0, nhcell-1 {
// HC -> GC
for j=0, 159 { // NOTE number of connections explicitly coded here!
Gauz3 = rdhc2gc.repick() // [-130:130]
//print Gauz3
if (i*83+41+Gauz3 > 499) {npost = i*83+41+Gauz3-500 } // NOTE: 500 is explicitly coded here!
if (i*83+41+Gauz3 < 0) {npost = i*83+41+Gauz3+500}
if ((i*83+41+Gauz3 > -1) && (i*83+41+Gauz3 < 500)) {npost = i*83+41+Gauz3}
//print npost
dbr = rdsyna.repick()
if (nhcell != 6) {max_conn = 5} else {max_conn = 3}
c1 = cells.object(npost)
c2 = cells.object(i+ngcell+nbcell+nmcell)
if ((is_connected(c1, c2) == 0) && (vhc2gc.x[npost] < max_conn)) {
nc_append(i+ngcell+nbcell+nmcell, npost, dbr+4, 0.5e-3*scale_HC2GC_strength, 1.6, 10)
vhc2gc.x[npost] +=1
if (debug_ == 1 ) {print i, npost, dbr+4}
} else {
j -= 1
if (debug_ == 1 ) {print "HC2GC"}
}
}
// HC -> BC
for j=0, 3 {
Gauz3 = rdhc2bc.repick() // [-2:2]
if (i+Gauz3 > nbcell-1) {npost = i+Gauz3-nbcell}
if (i+Gauz3 < 0) {npost = i+Gauz3+nbcell}
if ((i+Gauz3 > -1) && (i+Gauz3 < nbcell)) {npost = i+Gauz3}
dbr = rdsyna.repick()
c1 = cells.object(npost+ngcell)
c2 = cells.object(i+ngcell+nbcell+nmcell)
if ((is_connected(c1, c2) == 0) && (vhc2bc.x[npost] < 5)) { // NOTE < 5 coded explicitly here!
nc_append(i+ngcell+nbcell+nmcell, npost+ngcell, dbr+10, 0.5e-3, 1.6, 10) // Hcell to Bcell
if (debug_ == 1 ) {print npost, dbr+10}
vhc2bc.x[npost] += 1
} else {
j -= 1
if (debug_ == 1 ) {print "HC2BC"}
}
}
// HC -> MC
for j=0, 3 {
Gauz3 = rdhc2mc.repick() // [-2:2]
//print Gauz3
if (i*2+2+Gauz3 > 14) {npost = i*2+2+Gauz3-15 }
if (i*2+2+Gauz3 < 0) {npost = i*2+2+Gauz3+15}
if ((i*2+2+Gauz3 >-1) && (i*2+2+Gauz3 < 15)) {npost = i*2+2+Gauz3}
//print npost
dbr = rdsynb.repick()
// if ((is_connected(MossyCell[npost], HIPPCell[i]) == 0) && (vhc2mc.x[npost] < 2) && (killMC.contains(ngcell+nbcell+npost) == 0)) //use if killing MC
if (nhcell != 6) {max_conn = 4} else {max_conn = 2}
c1 = cells.object(npost+ngcell+nbcell)
c2 = cells.object(i+ngcell+nbcell+nmcell)
if ((is_connected(c1, c2) == 0) && (vhc2mc.x[npost] < max_conn)) { // NOTE < 2 coded explicitly here!! => changed to 4
nc_append(i+ngcell+nbcell+nmcell, npost+ngcell+nbcell, dbr+13, 1.5e-3, 1, 10) // Hcell to Mcell
//if (debug_ == 1 ) {print npost, dbr+13}
vhc2mc.x[npost] += 1
} else {
j -= 1
if (debug_ == 1 ) {print "HC2MC"}
}
// if (killMC.contains(ngcell+nbcell+npost) == 1){ j += 1 print "dead MC"} //use if killing MC
}
}
}
// ########################################################################
// FUNCTIONS FOR DATA STORAGE AND PRINTING
// ########################################################################
proc update_voltage_trace_matrix() { local i
/* Procedure called in every time step to add voltage of
recorded cell at time t
*/
VmT.append(t)
for i=0, (ngcell+nbcell +nmcell +nhcell -1) {
VmMat[i].append( cells.object[i].soma.v(0.5))
}
}
proc convert_voltage_trace_to_spikes() { local i, j, max_id
/* Converts membrane potential vector to a spike train
*/
max_id = npp
if ( print_template == 1) { max_id = npp }
if ( print_Raster == 1) { max_id = ngcell+nbcell +nmcell +nhcell-1 }
for i=0, (max_id) {
// Used to make a binary version of a spike train.
// Vector().spikebin(<data>, <thresh>)
// <data> is a vector of membrane potential.
// <thresh> is the voltage threshold for spike detection.
// <v> is set to all zeros except at the onset of spikes
// (the first dt which the spike crosses threshold)
Spike[i].spikebin(VmMat[i], 0)
}
}
obfunc myprecell() { localobj nil
/* Read out the ID of pre or precell (by Ted Carnevale)
Arguments:
$o1: NetCon
*/
// nil points to NULLobject
if ($o1.precell() == nil) {
return $o1.pre()
} else {
return $o1.precell()
}
}
proc save_network_connections(){ local i localobj dfile
/* Writes network adjacency list to file
*/
strdef DGNC_name_
dfile = new File()
sprint (DGNC_name_, "%s%s-%d-%d-%s.%s", datadir, "DGNC", PP_box_startid_, PP_box_nr_, idname, suffix)
print "\tWriting network connections to file: ", DGNC_name_
dfile.wopen(DGNC_name_)
for i=0, nclist.count-1 {
dfile.printf("%s\t%s\n", myprecell(nclist.object(i)), nclist.object(i).postcell)
}
dfile.close()
}
proc initialize_voltage_trace_file(){
/* Creates the header for the membrane voltage file
*/
strdef DGVt_name_ // Voltage file name
efile = new File()
// [ TODO ] -- n_to_print here is not good. Needs to be an object parameter or something.
n_to_print = 500
sprint (DGVt_name_,"%s%s-%d-%d-%s.%s", datadir, "DGVt", PP_box_startid_, PP_box_nr_, idname, suffix)
print "\tWriting membrane voltage traces of all cells to file: ", DGVt_name_
efile.wopen(DGVt_name_)
efile.printf("# t\t")
efile.printf("\n")
efile.close(DGVt_name_)
}
proc print_voltage_traces_to_file(){ local j
/* Prints voltage traces to file
TODO:
- Have better way to implement n_to_print so that
we can print not only granule cells, but still do this neatly
*/
efile.aopen(DGVt_name_)
efile.printf("%f\t", t)
for i = 0, n_to_print-1 {
efile.printf("%f\t", cells.object[int(i*ngcell/n_to_print)].soma.v(0.5))}
efile.printf("\n")
efile.close(DGVt_name_)
}
proc write_spike_times_to_file() { local i, j localobj spikefile
strdef DGsp_name_
sprint (DGsp_name_,"%s%s-%d-%d-%s.%s", datadir, "DGsp", PP_box_startid_, PP_box_nr_, idname, suffix)
print "\tWriting spike times to file: ", DGsp_name_
spikefile = new File()
spikefile.wopen(DGsp_name_)
k = 0
while (k < VmT.size) {
for j = 0, (ngcell+nbcell+nmcell+nhcell-1) {
if(Spike[j].x[k] != 0) {
// Writes out time of spike and cell id!
spikefile.printf("%f\t%d\n", VmT.x[k], j)
}
}
k +=1
}
spikefile.close(DGsp_name_)
}
proc write_stimuli_to_file() {local i, j localobj stimfile, noisestimfile
/* Writes stimuli to file
*/
stimfile = new File()
strdef f_name_
sprint (f_name_,"%s%s-%d-%d-%s.%s", datadir, "StimIn", PP_box_startid_, PP_box_nr_, idname, suffix)
print "\tWriting input trains to file: ", f_name_
stimfile.wopen(f_name_)
for i=0,npp-1 {
for j=0,vec_stim[i].size()-1 {
stimfile.printf("%f\t%d\n",vec_stim[i].x[j],i)
}
}
stimfile.close(f_name_)
// write out noise input...
if ((print_stim_noise == 1) && (spontaneous_activity_ == 1)) {
noisestimfile = new File()
strdef f_name_noise_
sprint (f_name_noise_,"%s%s-%d-%d-%s.%s", datadir, "NoiseStimIn", PP_box_startid_, PP_box_nr_, idname, suffix)
noisestimfile.wopen(f_name_noise_)
for i=0,ngcell-1 {
for j=0,vec_stim_noise[i].size()-1 {
noisestimfile.printf("%f\t%d\n",vec_stim_noise[i].x[j],i)
}
}
noisestimfile.close(f_name_noise_)
}
}
// ########################################################################
// FUNCTIONS TO INITIALIZE AND RUN SIMULATION
// ########################################################################
proc initialize_simulation_run() { local dtsav, temp, secsav
/* Initialize a simulation run
Arguments:
- $1 : int : 100 [ TODO ] NEED TO UNDERSTAND THIS MORE. SOMETHING ABOUT SEED SETTING
*/
finitialize(v_init)
// [ TODO ] - Figure out what these three next lines are for
t = -1000 // negative t step to initialize to steady state
dtsav = dt // original simulation step size
dt= 10 // larger step size
// if cvode is on, turn it off to do large fixed step
// [ TODO ] - Learn why
temp = cvode.active()
if ( temp != 0 ) { cvode.active(0) }
while( t < -100 ) { fadvance() }
//restore cvode if reqd
if (temp!=0) { cvode.active(1) }
dt = dtsav
t = 0
if ( cvode.active() ){
cvode.re_init()
} else {
fcurrent()
}
// restart number generators for PPStim
t = 0
finitialize(v_init)
trial_old = trial
if ($1 == 100) {
// Reset all Poisson inputs with the correct seed
for i = 0, rslist.count()-1 rslist.o(i).start(trial)
} else {
// reset int((rslist.count()-1)*(1-$1/100.)) Poisson inputs with a different seed
trial = trial + 1
for i = 0, int((rslist.count()-1)*(1-$1/100.)) rslist.o(i).start(trial)
// reset rslist.count()-int((rslist.count()-1)*(1-$1/100.)) Poisson inputs with the correct seed
trial = trial_old
for i = int((rslist.count()-1)*(1-$1/100.)), rslist.count()-1 rslist.o(i).start(trial)
}
// init noise Generators (with seed trial_noise ...)
// reset all Poisson noise inputs with the trial_noise seed
trial = trial_noise
for i = 0, rs_noiselist.count()-1 rs_noiselist.o(i).start(trial)
trial = trial_old
VmT = new Vector()
for i=0, cells.count-1 {
VmMat[i] = new Vector()
}
for i=0, (ngcell + nbcell + nmcell + nhcell -1) {
Spike[i] = new Vector() // vector of spikes for each cell
Spike_times[i] = new Vector() // vector of spikes for each cell
}
// finalize at the end
// [ TODO ] - Figure out what's going on here
t = 0
finitialize(v_init)
frecord_init()
// Select input pattern
for i=PP_box_startid_, PP_box_startid_+PP_box_nr_-1 {
PPSt[i].pp.start = PP_box_start_
PPSt[i].pp.forcestop = PP_box_stop_
if (PP_freq_) {
PPSt[i].pp.freq = PP_freq_
PPSt[i].pp.number = PP_box_nspk_
PPSt[i].pp.min_invl = PP_min_invl_
PPSt[i].pp.scale_max_invl = PP_scale_max_invl_
//print "pattern using osc"
} else {
PPSt[i].pp.nspk = PP_box_nspk_
PPSt[i].pp.status = 1
//print "pattern using box"
}
}
// Initialize the GC spontaneous activity generators
// [ TODO ] Need more general way to set the number that receive SA generators
// This needs to be based on the desired frequency of spont. activity.
if (spontaneous_activity_ == 1) {
// Rdsa_thresh is the number of granule cells that should be active
// with a single spike in the simulation window in order to achieve
// a spontaneous firing rate of `spontaneous_activity_rate`
rdsa_thresh = ngcell*(spontaneous_activity_rate/(1000/tstop))
for i=0,ngcell-1 {
rdsa_rv = rdsa.repick()
if (rdsa_rv <= rdsa_thresh) {
PPSt_noise[i].pp.status = 1
PPSt_noise[i].pp.start = 0.
PPSt_noise[i].pp.forcestop = tstop
PPSt_noise[i].pp.nspk = 1
}
if (rdsa_rv > rdsa_thresh) {
PPSt_noise[i].pp.status = 0
PPSt_noise[i].pp.start = 0.
PPSt_noise[i].pp.forcestop = tstop
PPSt_noise[i].pp.nspk = 0
}
}
}
}
proc run() {local i, j
// MAKE NETWORK
make_cells()
make_input_stimulation()
make_connections()
// INITIALIZE OPTIMIZER
cvode = new CVode()
// RE-INITIALIZE DATA CONTAINERS
objref Spike[ntotal-1], Spike_times[ntotal-1]
objref VmT, VmMat[cells.count]
// RUN SIMULATION
print "Running simulation."
// Write connections to file
save_network_connections()
for (PP_box_startid_ = 0; PP_box_startid_ <= n_patterns; PP_box_startid_ += 1 ) {
initialize_simulation_run(100)
// File header voltage output file
if (print_Vtrace == 1) {
initialize_voltage_trace_file()
}
while (t<tstop) {
fadvance()
if (print_Vtrace == 1) {
print_voltage_traces_to_file()
}
if ((print_Raster == 1) || (print_template == 1)) {
update_voltage_trace_matrix()
}
}
if ((print_Raster == 1) || (print_template == 1)) {
convert_voltage_trace_to_spikes()
}
if (print_Raster == 1) {
write_spike_times_to_file()
}
if (print_stim_==1) {
write_stimuli_to_file()
}
}
}
endtemplate DentateGyrus