Convergence regulates synchronization-dependent AP transfer in feedforward NNs (Sailamul et al 2017)

 Download zip file 
Help downloading and running models
Accession:241979
We study how synchronization-dependent spike transfer can be affected by the structure of convergent feedforward wiring. We implemented computer simulations of model neural networks: a source and a target layer connected with different types of convergent wiring rules. In the Gaussian-Gaussian (GG) model, both the connection probability and the strength are given as Gaussian distribution as a function of spatial distance. In the Uniform-Constant (UC) and Uniform-Exponential (UE) models, the connection probability density is a uniform constant within a certain range, but the connection strength is set as a constant value or an exponentially decaying function, respectively. Then we examined how the spike transfer function is modulated under these conditions, while static or synchronized input patterns were introduced to simulate different levels of feedforward spike synchronization. We observed that the synchronization-dependent modulation of the transfer function appeared noticeably different for each convergence condition. The modulation of the spike transfer function was largest in the UC model, and smallest in the UE model. Our analysis showed that this difference was induced by the different spike weight distributions that was generated from convergent synapses in each model. Our results suggest that the structure of the feedforward convergence is a crucial factor for correlation-dependent spike control, thus must be considered important to understand the mechanism of information transfer in the brain.
Reference:
1 . Sailamul P, Jang J, Paik SB (2017) Synaptic convergence regulates synchronization-dependent spike transfer in feedforward neural networks. J Comput Neurosci 43:189-202 [PubMed]
Citations  Citation Browser
Model Information (Click on a link to find other models with that property)
Model Type: Realistic Network; Synapse;
Brain Region(s)/Organism:
Cell Type(s): Hodgkin-Huxley neuron;
Channel(s): I Sodium; I Potassium; I T low threshold; I Cl, leak;
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: NEURON; MATLAB;
Model Concept(s): Synchronization; Oscillations; Action Potentials; Activity Patterns; Information transfer; Synaptic Convergence;
Implementer(s): Sailamul, Pachaya [pachaya_sailamul at brown.edu]; Jang, Jaeson ; Paik, Se-Bum ;
Search NeuronDB for information about:  I T low threshold; I Sodium; I Potassium; I Cl, leak;
////////////////////////////////////////////////////////////////////////////////////////////////////
////////// Convergence Connections in Feedforward network L1 = Poisson Generator Cells , L2 = E cells
////////////////////////////////////////////////////////////////////////////////////////////////////

load_file("nrngui.hoc") //load basics library
T_everythng = startsw()
load_file("ranstream.hoc") //for InGauss random seed
load_file("Parameters.hoc") //Parameters
load_file("CellsTemplate.hoc") //load basics library/cells' template  Note: Parameters for Cells template are specified inside

////////////////////////////////////////////////////////////////////////////////////////////////////
////////// Variables initialization
////////////////////////////////////////////////////////////////////////////////////////////////////
// These are the default values
// They may get overwrite when the simulation is called at the end of this script

// Constant
RANDOM123_ID1_POISSONSPK = 1 // ID 1 = feed forward input from L1 to L2

// Simulation Setting
//CELL_TYPE = 1 -----------------------------------------------------------------------
TRIAL_NO = 1 // Initialization purpose only
CONVERGENT_TYPE = 1 // Initialization purpose only

strdef CnvrgentConnTypeTxt 
proc get_CnvrgentConnTypeTxt() {

if(CONVERGENT_TYPE ==1){ CnvrgentConnTypeTxt  = "GG"} 
if(CONVERGENT_TYPE ==2){ CnvrgentConnTypeTxt  = "UC"} 
if(CONVERGENT_TYPE ==3){ CnvrgentConnTypeTxt  = "UE"} 

} 
get_CnvrgentConnTypeTxt()

// Connection
W_SCALE = 1e-05  // since the value of W is really small, a scaling factor is used,the "weightingFactor" is specified relative to this scale
weightingFactor = 50 // Weighting Factor or connection strength is weightingFactor x W_SCALE
range = 50 //Range of connection


// Pattern of Input
OSC_F =  40 // Oscillation frequency
PHASE = 0 // Input Phase
OSC_rltAmp = 0// [0, 0.5, 1] //Relative amplitude for each synchronization level of input
Input_spk_avg_fr = 20 // The average firing rate
L1_PHASE_RANDOM_SIG = 0 // 0, 0.25 0.5 0.75 1 --------> Level of heterogeneity in input pattern


//Initialized values
steps_per_ms = 1    // Resolution for data gathering set
v_init = -70 // initial membrane potential
tstop = 5000  // simulation runtime

// Simulation Code
strdef SimCode, SimCtrl
SimCtrl="ModelDB"
//Generate simulation code for current set of parameters

proc getSimCode(){
//sprint(SimCode, "%s_%s_InputFR%g_OscF%g_OscrltAmp%g_OscPhaseSig%g_Wscale%.5f_W%g_range%g_Trial%g_T%g", SimCtrl, CnvrgentConnTypeTxt, Input_spk_avg_fr, OSC_F, OSC_rltAmp, L1_PHASE_RANDOM_SIG, W_SCALE, weightingFactor, range,TRIAL_NO, tstop) //   
sprint(SimCode, "%s_%s_InputFR%g_OscF%g_OscrltAmp%g_Wscale%.5f_W%g_range%g_Trial%g_T%g", SimCtrl, CnvrgentConnTypeTxt, Input_spk_avg_fr, OSC_F, OSC_rltAmp, W_SCALE, weightingFactor, range,TRIAL_NO, tstop) //   
}
//getSimCode()

proc getSimCode_hg(){
sprint(SimCode, "%s_%s_InputFR%g_OscF%g_OscrltAmp%g_OscPhaseSig%g_Wscale%g_W%g_range%g_Trial%g_T%g", SimCtrl, CnvrgentConnTypeTxt, Input_spk_avg_fr, OSC_F, OSC_rltAmp, L1_PHASE_RANDOM_SIG, W_SCALE, weightingFactor, range,TRIAL_NO, tstop) //   

}
//getSimCode_hg()

// NOTE : How to reuse List and Vector  --> IClamplist.remove_all(),  p.resize(0)

////////////////////////////////////////////////////////////////////////////////////////////////////
//////////  Setting Directories for input and output
////////////////////////////////////////////////////////////////////////////////////////////////////

strdef dirResults, dirInFiles, dirPhaseFiles
dirResults ="SimResult/" //Directory for Simulation Results
dirInFiles = "Input/" //Location of feedforward connection
dirPhaseFiles = "Heterogeneity/"
////////////////////////////////////////////////////////////////////////////////////////////////////
////////// Network specification interface (Helper function)
////////////////////////////////////////////////////////////////////////////////////////////////////

objref cells, cellsIN, cellsE, cellsI, cellsCN, nclist, netcon, cellsFFin 
{cells = new List() cellsIN = new List() cellsE = new List() cellsI = new List()  cellsCN = new List()
 nclist = new List() cellsFFin= new List() }

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

func nc_append() { local w, delay 
//srcindex, tarcelindex, synindex, weight, delay   //Ex. //  /* C1 -> C2.E0 */  nc_append(1,   2, 0,  0.12,1)
  if ($3 >= 0) {
    netcon = cells.object($1).connect2target(cells.object($2).synlist.object($3))  // Excitatory and Inhibitory effects define at the synaptic input of target cell (*post* synaptic)
    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) //nclist is list of NetCon object, The source cell can be access by call netcon.precell and access target cell by calling netcon.postcell ex. nclist.o(1).precell 
  return nclist.count - 1 
}
printf("Done Network specification interface")

////////////////////////////////////////////////////////////////////////////////////////////////////
////////// Download Cells Position
////////////////////////////////////////////////////////////////////////////////////////////////////
strdef layer1locTxt, layer2locTxt 
layer1locTxt = "Cells_position_source_layer1.txt" // Source Layer
layer2locTxt = "Cells_position_target_layer2.txt"// Target Layer

//load position files : Layer 1 
objref posL1x, posL1y, posL2x, posL2y
objref fin
fin = new File()
fin.ropen(layer1locTxt)  
 
nL1 = fin.scanvar()
posL1x = new Vector(nL1)
posL1y = new Vector(nL1)

for (i=0;i<nL1;i=i+1){
	posL1x.x(i) = fin.scanvar()
	posL1y.x(i) = fin.scanvar()
}
fin.close() 

//load position files : Layer 2
fin = new File()
fin.ropen(layer2locTxt)  
 
nL2 = fin.scanvar()
posL2x = new Vector(nL2)
posL2y = new Vector(nL2)

for (i=0;i<nL2;i=i+1){
	posL2x.x(i) = fin.scanvar()
	posL2y.x(i) = fin.scanvar()
}
fin.close() 

////////////////////////////////////////////////////////////////////////////////////////////////////
////////// Create Cells
////////////////////////////////////////////////////////////////////////////////////////////////////
//// Create the cells in Layer 2 (H-H model)
objref cellsL2
cellsL2 =  new List()

cIDrec = -1 //cell counter
	str_E = cIDrec+1  
	L2rec= -1 //cell counter
	for (i=0;i<nL2;i=i+1){
		cIDrec = cIDrec+1 
		L2rec=L2rec+1
		cell_append(new Target_Cell(),posL2x.x(i),posL2y.x(i), 0,cIDrec, L2rec,1)
		cellsL2.append(cells.object(cIDrec))
	} //End generate L2

//// Create the cell in Layer 1 , for Poisson Spike Generator 
objref cellsL1
cellsL1 = new List()

	str_FFin = cIDrec+1
	FFinrec = -1 //Cell counter
	nnInspkCell = nL1
	FFzpos = -50 // Position of cell in z-axis - Our cells are located in x-y plane, thus we set position in z-axis to an arbitrary value
	for (i=0;i < nnInspkCell; i=i+1){
		cIDrec = cIDrec+1 
		FFinrec=FFinrec+1
		cell_append(new In_spk_VecStim(),posL1x.x(i),posL1y.x(i), FFzpos,cIDrec, FFinrec,0) // cell type = 0  for FF input
		cellsL1.append(cells.object(cIDrec))
	}
	
////////////////////////////////////////////////////////////////////////////////////////////////////
////////// Specified Function for Poisson Generator (Input Spike train)
////////////////////////////////////////////////////////////////////////////////////////////////////

obfunc  OscilatingInputFR() {local osc_f, phase, osc_amp, meanFR,time localobj v1
// FR = A sin(2*pi*f + phase) ---- equation
// $1 = osc_f, $2 = phase (in rad), $3 = osc_amp, $4 = meanFR, $5 = resolution/ steps_per_ms , $6 = time = tstop 
{ osc_f =$1  phase =$2 osc_amp=$3 meanFR=$4 }

time = $6/ $5
v1 = new Vector(time)
v1.sin( osc_f,phase, 1/ $5) // dt = 1 ms
v1.mul(osc_amp)
v1.add(meanFR) 
return v1 // v1 = expected firing rate at each time bin
} 

obfunc  poissonGenerator() { local tmpr,ii localobj rr, spktrain, spktime, p
//Return vector of spiking time (in ms)
// $1 = seed, $2 = resolution (size of one bin in ms) = steps_per_ms, $3 = time (ms), $4 = average firing rate
	rr = new Random()
	rr.uniform(0,1)
	rr.Random123(RANDOM123_ID1_POISSONSPK,$1,TRIAL_NO)
	spktrain = new Vector($3/$2) 
	spktime = new Vector()	
	p = OscilatingInputFR(OSC_F, PHASE,$4*OSC_rltAmp,$4, $2, $3) // // $1 = osc_f, $2 = phase (in rad), $3 = osc_amp, $4 = meanFR, $5 = resolution/ steps_per_ms , $6 = time = tstop 
	p.div(1000/$2) // chance of spike to happen (1 ms resolution)
	
	for ii =0,spktrain.size-1 {
		tmpr = rr.repick
		
		if (tmpr < p.x(ii)){  //spike occur
			spktrain.x(ii) = 1
			spktime.append($2*ii)
		}else{
			spktrain.x(ii) = 0
		}
	}
	p.resize(0)
	return spktime // spkvec, 1) the multiple arrival time in one VecStim is account as the only one spike. 2) The spike time vector need to be sort ascending(less...more).
}
 
//poissonGenerator_hg = poissonGenerator with heterogenity in each cell (variety in phase)
obfunc  poissonGenerator_hg() { local tmpr,ii localobj rr, spktrain, spktime, p
//Return vector of spiking time (in ms)
// $1 = seed, $2 = resolution (size of one bin in ms) = steps_per_ms, $3 = time (ms), $4 = average firing rate, $5 = phase
	rr = new Random()
	rr.uniform(0,1)
	rr.Random123(RANDOM123_ID1_POISSONSPK,$1,TRIAL_NO)
	spktrain = new Vector($3/$2) 
	spktime = new Vector()	
	p = OscilatingInputFR(OSC_F, $5,$4*OSC_rltAmp,$4, $2, $3) // // $1 = osc_f, $2 = phase (in rad), $3 = osc_amp, $4 = meanFR, $5 = resolution/ steps_per_ms , $6 = time = tstop 
	p.div(1000/$2) // chance of spike to happen (1 ms resolution)
	
	for ii =0,spktrain.size-1 {
		tmpr = rr.repick
		
		if (tmpr < p.x(ii)){  //spike occur
			spktrain.x(ii) = 1
			spktime.append($2*ii)
		}else{
			spktrain.x(ii) = 0
		}
	}
	p.resize(0)
	return spktime // spkvec, 1) the multiple arrival time in one VecStim is account as the only one spike. 2) The spike time vector need to be sort ascending(less...more).
}
 
////////////////////////////////////////////////////////////////////////////////////////////////////
////////// Add Poisson Spike Generator to L1
////////////////////////////////////////////////////////////////////////////////////////////////////
objref spkvec_list
spkvec_list = new List()

objref recordseed
recordseed = new Vector() 
recordseed.resize(cellsL1.count)

proc AddInputSpikeVector_to_L1() {
spkvec_list.remove_all()

recordseed = new Vector() 
recordseed.resize(cellsL1.count)

for i=0, cellsL1.count-1 { // Record Poisson Spike Train for each L1 cell

spkvec_list.append( poissonGenerator(cellsL1.o(i).cID,steps_per_ms,tstop,Input_spk_avg_fr ))
recordseed.x(i) = cellsL1.o(i).cID // specified ID#2 of RANDOM123 as cell ID of that input cell--> to be sure that there is no any two cells with same seed
cellsL1.o(i).pp.play(spkvec_list.o(i)) 
	
}
} 
//AddInputSpikeVector_to_L1()
strdef PHASEfile
proc AddInputSpikeVector_to_L1_hg() {local spkPhase, nnPhaseL1 localobj fin
spkvec_list.remove_all()

recordseed = new Vector() 
recordseed.resize(cellsL1.count)
fin = new File()

//// For each heterogeneity level(sigRand)
// Read phase list file : Heterogeneity_N1150_RandomSig0.5_Trial1
sprint(PHASEfile, "%sHeterogeneity_N%g_RandomSig%g_Trial%g.txt", dirPhaseFiles,nL1, L1_PHASE_RANDOM_SIG, TRIAL_NO)    
//Ex : Heterogeneity_N1150_RandomSig0.5_Trial1
fin.ropen(PHASEfile) 
nnPhaseL1  = fin.scanvar()

////////////////////////////////////////////

for i=0, cellsL1.count-1 { // Record Poisson Spike Train for each L1 cell
spkPhase = fin.scanvar()
spkvec_list.append( poissonGenerator_hg(cellsL1.o(i).cID,steps_per_ms,tstop,Input_spk_avg_fr,spkPhase ))
recordseed.x(i) = cellsL1.o(i).cID // specified ID#2 of RANDOM123 as cell ID of that input cell--> to be sure that there is no cell with same seed
cellsL1.o(i).pp.play(spkvec_list.o(i)) 
	
}

fin.close(PHASEfile)  
} 
//AddInputSpikeVector_to_L1_hg()

////////////////////////////////////////////////////////////////////////////////////////////////////
//////////  Generate L1 - L2 connection
////////////////////////////////////////////////////////////////////////////////////////////////////
strdef  FRfile
objref CorrInspk_info,srcV, tarV, wVec, dVec
{CorrInspk_info = new Vector() srcV = new Vector() tarV = new Vector() wVec = new Vector()
dVec = new Vector()}
//strdef dirInFiles
//dirInFiles = "Input/"
proc MakeL1_L2_Conn(){local nnConnL1,nnInspkCell localobj fin, fout
nclist.remove_all()
{CorrInspk_info.resize(0) srcV.resize(0)	tarV.resize(0)	wVec.resize(0)	dVec.resize(0)}
fin = new File()

// Read Connection Files  
sprint(FRfile, "%sConvergentInput_%s_Wscale%.5f_W%g_range%g_Trial%g.txt", dirInFiles,CnvrgentConnTypeTxt, W_SCALE, weightingFactor, range,TRIAL_NO)    
//Ex : ConvergentInput_GG_Wscale0.00001_W50_range50_Trial1
fin.ropen(FRfile) 
nnConnL1  = fin.scanvar()
nnInspkCell  = fin.scanvar() 


CorrInspk_info = new Vector(nnConnL1)

srcV = new Vector(nnConnL1)
tarV = new Vector(nnConnL1)
wVec = new Vector(nnConnL1)
dVec = new Vector(nnConnL1)

	for (i=0;i<nnConnL1 ;i=i+1){
		srcV.x(i) = fin.scanvar()
		tarV.x(i) = fin.scanvar()
		wVec.x(i) = fin.scanvar()
		dVec.x(i) = fin.scanvar()
		
		// Make Connection				
		nc_append(cellsL1.o(srcV.x(i)).cID,cellsL2.o(tarV.x(i)).cID,0,wVec.x(i), dVec.x(i)) //srcindex, tarcelindex, synindex, weight, delay 
	}
fin.close(FRfile)
	
fout = new File()
sprint(FRfile, "%sRecordConnList_%s.txt", dirResults ,SimCode)
fout.wopen(FRfile) 


fout.printf("%g\t%g\t%g\t%g\n",nL1, nL2, nclist.count, tstop) //No of cell in layer 1 , No of cell in layer 2, number of all connections
for i = 0, nclist.count-1{
	fout.printf("%g\t%g\t%g\t%g\n",srcV.x(i), tarV.x(i), wVec.x(i), dVec.x(i))	 // presynaptic cell ID, postsynaptic cell ID, connection strength, delay(= time at which spike from source cell arrive at target cell)
}
fout.close()
} 
//MakeL1_L2_Conn()
	
////////////////////////////////////////////////////////////////////////////////////////////////////
//////////  Make Vector to record L2 cell activity 
////////////////////////////////////////////////////////////////////////////////////////////////////

objref vVec_List, i_cap_List
vVec_List= new List() 
i_cap_List = new List() 

objref vVec, i_capVec
vVec= new Vector() 
i_capVec= new Vector() 

vVec_List.remove_all() 
i_cap_List.remove_all() 


//Record all V in all cells 
for id=0, cellsL2.count-1 {
	vVec = new Vector()
	i_capVec = new Vector()
	
	vVec.record(&cellsL2.o(id).soma.v(0.5),RESOLUTION) //record soma's voltage with the resolution of 1 ms
	i_capVec.record(&cellsL2.o(id).soma.i_cap(0.5),RESOLUTION)
	
	vVec_List.append(vVec)
	i_cap_List.append(i_capVec)	
}


////////////////////////////////////////////////////////////////////////////////////////////////////
//////////  Run Simulation
////////////////////////////////////////////////////////////////////////////////////////////////////

proc run_all() {
trun = startsw()
finitialize()
frecord_init()
print "Called run_all()"
tstop = $1  
steps_per_ms = 1   
v_init = -70
run()
print "Finished run_all()"
print "Total Run Time ", startsw() - trun
}
// run_all(tstop)

////////////////////////////////////////////////////////////////////////////////////////////////////
//////////  Save Neuron Activity to files
////////////////////////////////////////////////////////////////////////////////////////////////////
// save L1: spikes train for each cell, L2 : membrane potential and membrane current

strdef saveVecFName
saveVecFName = ""

proc save_vectors_to_file(){local i,j localobj fout,id
fout = new File()
fout.wopen(saveVecFName)
fout.printf("%g\t%g\t%g",cellsL2.count,0, tstop) //#E, #I, tstop[ms]  (E cells' ID always go first) :: there is no inhibitory neuron(I) in this network
for i = 0,$o1.count-4{
	fout.printf("\t%g",0)
}
fout.printf("\n")

id = new Vector()
id.indgen(0, $o1.count-1,1)
id.printf(fout,"%g\t",0,id.size-1)
//fout.printf("%g\n",id.x(id.size-1))

for i = 0, $o1.o(0).size-1{
	for j = 0, $o1.count-2 {
		fout.printf("%g\t",$o1.o(j).x(i))
	}
	fout.printf("%g\n",$o1.o($o1.count-1).x(i))
}
fout.close()
}


proc save_SPKtrain_to_file(){local i,j localobj fout,id

fout = new File()
fout.wopen(saveVecFName)
fout.printf("%g\t%g\t%g\n",cellsL1.count,0, tstop) //#E, #I, tstop[ms]  

for i = 0, spkvec_list.count-1{
	fout.printf("%g\t",i )
	if(spkvec_list.o(i).size >0){
	spkvec_list.o(i).printf(fout,"%g\t",0, spkvec_list.o(i).size-1)
	}
	fout.printf("\n")
}
fout.close()

}


// save file  
proc save_all(){
trun = startsw()

//L1 Spikes train
sprint(saveVecFName,"%sInputSpkTrain_%s.txt",dirResults,SimCode) 
save_SPKtrain_to_file() // Input spike train

//L2 Voltage
sprint(saveVecFName,"%sSomaVolt_%s.txt",dirResults,SimCode)
save_vectors_to_file(vVec_List) // membrane potential of target cells at soma

//L2 i_cap
sprint(saveVecFName,"%sSoma_i_cap_%s.txt",dirResults,SimCode) 
save_vectors_to_file(i_cap_List) //capacitance of target cells
print "Total Time for Saving vectors to files : ", startsw() - trun
}
//save_all()


////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////// Simulations
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
proc new_Loop_Run(){
//print " ==================================================== "
print "SIMULATION CODE = ", SimCode
//print " ==================================================== "
AddInputSpikeVector_to_L1()
MakeL1_L2_Conn()
run_all(tstop)
save_all()
}

proc new_Loop_Run_hg(){
//print " ==================================================== "
print "SIMULATION CODE = ", SimCode
//print " ==================================================== "
AddInputSpikeVector_to_L1_hg()
MakeL1_L2_Conn()
run_all(tstop)
save_all()
}
////////////////////////////////////////////////////
// Network Configuration
////////////////////////////////////////////////////
/*
Parameters List
    CnvrgentConnTypeTxt : Code for convergence connection model (One of 'GG', 'UC', 'UE')
    Input_spk_avg_fr: Input average firing rate (Default: 20Hz)
    OSC_F: Oscillation Frequency of source layer (Default: 40Hz)
    OSC_rltAmp: Oscillation relative amplitude (Strong = 1, Weak = 0.5, static = 0)
    weightingFactor: weight for connection strength
    range: range of connection
    TRIAL_NO: Trial number
    tstop : Simulation time (Default: 5000)
Return
*/

// Directories
dirInFiles = "Input/" //Location of feedforward connections information
dirPhaseFiles = "Heterogeneity/"
dirResults ="SimResult/" //Directory for Simulation Results

objref RANGE_LST , W_LST, OSC_AMP_LST, L1_PHASE_RANDOM_SIG_LST // Range and weighting factor are the two parameters that control level of convergence
RANGE_LST  = new Vector(1) // RANGE 
RANGE_LST.x(0) = 80 
range = RANGE_LST.x(0)

W_LST = new Vector(1) // WEIGHTING FACTOR
W_LST.x(0) =80
weightingFactor = W_LST.x(0)

OSC_AMP_LST = new Vector() // OSCILLATION AMPLITUDE (Strong = 1, Weak = 0.5, static = 0)
OSC_AMP_LST.indgen(0,1,0.5)  //0 0.5 1

L1_PHASE_RANDOM_SIG_LST = new Vector()
L1_PHASE_RANDOM_SIG_LST.indgen(0,1,0.25)

HETEROGENEITY_TEST = 0

// Note : tstop have to specified in above "Parameters" part under "//Initialized values" (Line#51 )
NUM_CONVERGENT_TYPE = 3 // Convergent Type 1:GG , 2: UC, 3: UE
N_Trial = 10
tstop = 100
SimCtrl="ModelDB"

// Without Heterogeneity testing

for TRIAL_NO = 1, N_Trial {	
	printf("==========================================================\n")
	printf("\tTRIAL NO = %g\n" , TRIAL_NO)
	printf("==========================================================\n")
	
    for o_ii =0, OSC_AMP_LST.size-1 {
		OSC_rltAmp = OSC_AMP_LST.x(o_ii)
		//printf("\t==========================================================\n")
		printf("\t\t Oscillation relative amplitude = %g\n" , OSC_rltAmp)
		printf("\t==========================================================\n")

		for t_ii = 1, NUM_CONVERGENT_TYPE { // Innermost
			CONVERGENT_TYPE = t_ii
			get_CnvrgentConnTypeTxt()
			getSimCode()
			printf("\t\t\t\t Convergent Connection Type  = %s\n" , CnvrgentConnTypeTxt)
			printf("\t\t\t==========================================================\n")
			new_Loop_Run()
			}
	}
}

// With Heterogeneity testing
if(HETEROGENEITY_TEST){
    for TRIAL_NO = 1, N_Trial {	
        printf("==========================================================\n")
        printf("\tTRIAL NO = %g\n" , TRIAL_NO)
        printf("==========================================================\n")
        
        for phi_ii = 0 , L1_PHASE_RANDOM_SIG_LST.size-1 {
            L1_PHASE_RANDOM_SIG = L1_PHASE_RANDOM_SIG_LST.x(phi_ii)
            //printf("\t==========================================================\n")
            printf("\t\t L1_PHASE_RANDOM_SIG = %g\n" , L1_PHASE_RANDOM_SIG)
            printf("\t==========================================================\n")
            
            for o_ii =0, OSC_AMP_LST.size-1 {
                OSC_rltAmp = OSC_AMP_LST.x(o_ii)
                //printf("\t==========================================================\n")
                printf("\t\t Oscillation relative amplitude = %g\n" , OSC_rltAmp)
                printf("\t==========================================================\n")
                
                for t_ii = 1, NUM_CONVERGENT_TYPE { // Innermost
                    CONVERGENT_TYPE = t_ii
                    get_CnvrgentConnTypeTxt()
                    getSimCode_hg()
                    printf("\t\t\t\t Convergent Connection Type  = %s\n" , CnvrgentConnTypeTxt)
                    printf("\t\t\t==========================================================\n")
                    new_Loop_Run_hg()
                }
            }
        }
    }

}


////////////////////////////////////////////////////////////////////////////////////////////////////
////////// End
////////////////////////////////////////////////////////////////////////////////////////////////////

print "Time for everything = ",  startsw() - T_everythng