// Simulate background activity 
// Papoutsi, 28/08/2014
// Modified by Petousakis K-E. to allow for the removal (i.e. lack of allocation) of synapses, 2019

//Vectors of Spike trains
objref BG_basal[syn_basal+2], BG_apical[syn_apical+2], BG_Inh_soma[syn_inh_soma+2], BG_Inh_apical[syn_inh_apical+2], BG_Inh_basal[syn_inh_basal+2]
//Background synapses
objref bc_apical_ampa[syn_apical+2], bc_apical_nmda[syn_apical+2], bc_basal_ampa[syn_basal+2], bc_basal_nmda[syn_basal+2], bc_gaba_soma[syn_inh_soma+2], bc_gaba_apical[syn_inh_apical+2], bc_gaba_basal[syn_inh_basal+2]
//VecStims
objref vc_apical[syn_apical+2], vc_basal[syn_basal+2], vc_gaba_soma[syn_inh_soma+2], vc_gaba_apical[syn_inh_apical+2], vc_gaba_basal[syn_inh_basal+2]
//NetCons
objref nc_apical_ampa[syn_apical+2], nc_apical_nmda[syn_apical+2], nc_basal_ampa[syn_basal+2], nc_basal_nmda[syn_basal+2], nc_gaba_soma[syn_inh_soma+2], nc_gaba_apical[syn_inh_apical+2], nc_gaba_basal[syn_inh_basal+2]
//Random objects
objref r, r_time, r_time_inh

//Make distribuiton of synapses according to length
objref mat_bas, mat_ap, mat_inh_bas, mat_inh_ap,vd
mat_bas=new Vector()
mat_ap=new Vector()
mat_inh_bas=new Vector()
mat_inh_ap=new Vector()
tot_ap=0
tot_bas=0
for i=0, apical.count()-1 {tot_ap=tot_ap+apical.o(i).sec.L}
for i=0, apical.count()-1 {
	mat_ap.append(int((apical.o(i).sec.L/tot_ap)*syn_apical))
	mat_inh_ap.append(int((apical.o(i).sec.L/tot_ap)*syn_inh_apical))
}
for i=0, basal.count()-1 {tot_bas=tot_bas+basal.o(i).sec.L}
for pd=0,4 {
	for i=0, list_primary[pd].count()-1 {
		mat_bas.append(round((list_primary[pd].o(i).sec.L/tot_bas)*syn_basal))
		mat_inh_bas.append(round((list_primary[pd].o(i).sec.L/tot_bas)*syn_inh_basal))
	}	
}


proc background_activity () {	
	//Define random for position
	r=new Random(num_seed2)
	r.uniform(0,1)
	//and for poisson trains
	r_time=new Random(num_seed+1200)
	r_time.poisson(Hz/1000)

	r_time_inh=new Random(num_seed-7000)
	r_time_inh.poisson(Hz_inh/1000)
	//Basal	excitatory background
	if (!nbB) {   //If nbB (no background basal) is 1, do not allocate these synapses
		syn=0
		c_d=-1
		for pd=0,4 {
			for num=0, list_primary[pd].count()-1 {
				c_d=c_d+1
				for i=0, mat_bas.x(c_d)-1 {
					BG_basal[syn]=new Vector()
					for tt=0, tstop-1 {
						if (r_time.repick()){
							BG_basal[syn].append(tt)
						}
					}
					vc_basal[syn] = new VecStim(0.5)
					vc_basal[syn].delay = 0
					vc_basal[syn].play(BG_basal[syn])
					PID=r.repick()
					list_primary[pd].o(num).sec bc_basal_ampa[syn]=new GLU(PID)
					list_primary[pd].o(num).sec bc_basal_nmda[syn]=new nmda(PID)
					if (!cut_basal) {
						// flagBG from L23V1_model_sim.hoc controls hard-coded nullification of basal synapses.
						// Value is 0, so the second "if" statement always runs. Rest are there for testing purposes.
						if(pd==0 && flagBG==1){  
							nc_basal_ampa[syn] = new NetCon(vc_basal[syn], bc_basal_ampa[syn], -20, 0, ampa_g*0)
							nc_basal_nmda[syn] = new NetCon(vc_basal[syn], bc_basal_nmda[syn], -20, 0, nmda_g*0)
						}
						if(pd!=0 || flagBG==0){
							nc_basal_ampa[syn] = new NetCon(vc_basal[syn], bc_basal_ampa[syn], -20, 0, ampa_g)
							nc_basal_nmda[syn] = new NetCon(vc_basal[syn], bc_basal_nmda[syn], -20, 0, nmda_g)
						}
						
					} else if (pd !=dend_id1 && pd!=dend_id2) {
						nc_basal_ampa[syn] = new NetCon(vc_basal[syn], bc_basal_ampa[syn], -20, 0, ampa_g)
						nc_basal_nmda[syn] = new NetCon(vc_basal[syn], bc_basal_nmda[syn], -20, 0, nmda_g)
					}
					syn=syn+1
				}
			}
		}
	}
	//Apical excitatory background
	if (!nbA) {   //If nbA (no background apical) is 1, do not allocate these synapses
		syn=-1
		for num=0, apical.count()-1 {
			for i=0, mat_ap.x(num)-1 {
				syn=syn+1
				BG_apical[syn]=new Vector()
				for tt=0, tstop-1 {
					if (r_time.repick()){
						BG_apical[syn].append(tt)
					}
				}

				vc_apical[syn] = new VecStim(0.5)
				vc_apical[syn].delay = 0
				vc_apical[syn].play(BG_apical[syn])

				PID=r.repick()
				apical.o(num).sec bc_apical_ampa[syn]=new GLU(PID)
				apical.o(num).sec bc_apical_nmda[syn]=new nmda(PID)

				if (!ablated) {
					nc_apical_ampa[syn] = new NetCon(vc_apical[syn], bc_apical_ampa[syn], -20, delFB_bg_exc, ampa_g)
					nc_apical_nmda[syn] = new NetCon(vc_apical[syn], bc_apical_nmda[syn], -20, delFB_bg_exc, nmda_g)
				}
			}
		}

	}

	// Inhibition Soma
	if (!niS) {

		for syn=0,syn_inh_soma-1 {
			BG_Inh_soma[syn]=new Vector()
			for tt=0, tstop-1 {
				if (r_time_inh.repick()){
					BG_Inh_soma[syn].append(tt)
				}
			}
			vc_gaba_soma[syn] = new VecStim(0.5)
			vc_gaba_soma[syn].delay = 0
			vc_gaba_soma[syn].play(BG_Inh_soma[syn])
			
			PID=r.repick()
			soma bc_gaba_soma[syn]=new GABAa(PID)
			nc_gaba_soma[syn] = new NetCon(vc_gaba_soma[syn], bc_gaba_soma[syn], -20, 0, gaba_g)
		}
	}


	//Inhibition Basal	
	if (!niB) {  //If niB (no inhibition basal) is 1, then we don't want this part to execute, or an error might occur.
		syn=0
		c_d=-1

		for pd=0,4 {  
			for num=0, list_primary[pd].count()-1 {
				c_d=c_d+1
				for i=0, mat_inh_bas.x(c_d)-1 {
					BG_Inh_basal[syn]=new Vector()
					for tt=0, tstop-1 {
						if (r_time_inh.repick()){
							BG_Inh_basal[syn].append(tt)
						}
					}
					vc_gaba_basal[syn] = new VecStim(0.5)
					vc_gaba_basal[syn].delay = 0
					vc_gaba_basal[syn].play(BG_Inh_basal[syn])
					
					PID=r.repick()

					list_primary[pd].o(num).sec bc_gaba_basal[syn]=new GABAa(PID)
					if (!cut_basal) {
						nc_gaba_basal[syn] = new NetCon(vc_gaba_basal[syn], bc_gaba_basal[syn], -20, 0, gaba_g)
					} else if (pd !=dend_id1 && pd!=dend_id2) {
						nc_gaba_basal[syn] = new NetCon(vc_gaba_basal[syn], bc_gaba_basal[syn], -20, 0, gaba_g)
					}
					syn=syn+1
				}
			}

		}
	}
	//Inhibition Apical
	if (!niA) {   //If niA (no inhibition apical) is 1, then we don't want this part to execute, or an error might occur.
		syn=-1
		for num=0, apical.count()-1  {	
			for i=0, mat_inh_ap.x(num)-1 {
				syn=syn+1
				BG_Inh_apical[syn]=new Vector()
				for tt=0, tstop-1 {
					if (r_time_inh.repick()){
						BG_Inh_apical[syn].append(tt)
					}
				}
				vc_gaba_apical[syn] = new VecStim(0.5)
				vc_gaba_apical[syn].delay = 0
				vc_gaba_apical[syn].play(BG_Inh_apical[syn])

				PID=r.repick()

				apical.o(num).sec bc_gaba_apical[syn]=new GABAa(PID)
				if (!ablated) {
					nc_gaba_apical[syn] = new NetCon(vc_gaba_apical[syn], bc_gaba_apical[syn], -20, delFB_bg_inh, gaba_g)
				}
			}
		}	
	}

// Commented-out code is for synaptic count validation, unnecessary for simulation execution.
/*	
	sbgexcB=0
	sbgexcA=0
	sbginhB=0
	sbginhA=0
	for i=0,6{
		sbgexcB+=mat_bas.x(i)
		sbginhB+=mat_inh_bas.x(i)
		print "Basal dendrite ",i, " receives ", mat_bas.x(i), " background excitatory synapses."
		print "Basal dendrite ",i, " receives ", mat_inh_bas.x(i), " background inhibitory synapses."
	}
	for i=0,42{
		sbgexcA+=mat_ap.x(i)
		sbginhA+=mat_inh_ap.x(i)
		print "Apical dendrite ",i, " receives ", mat_ap.x(i), " background excitatory synapses."
		print "Apical dendrite ",i, " receives ", mat_inh_ap.x(i), " background inhibitory synapses."
	}
	print "Actual Synaptic Totals: "
	print "Background Excitatory Basal: ", sbgexcB, " (defined as ", syn_basal, " )"
	print "Background Inhibitory Basal: ", sbginhB, " (defined as ", syn_inh_basal, " )"
	print "Background Excitatory Apical: ", sbgexcA, " (defined as ", syn_apical, " )"
	print "Background Inhibitory Apical: ", sbginhA, " (defined as ", syn_inh_apical, " )"
*/	
}