// written by Gabriela Cirtala
// Heterogenous ion channel conductance density model
// This code corresponds to Branch 7 of the PC (Figure 4G-I in Cirtala2024)
// This code only runs the case of 42PF (i.e. the threshold for br7)
// Shows the dendritic spike propagation (check the shape plot at the PF activation time)

// this code is based on the homogenous model proposed by Zang et al. (Zang 2021).

load_file("nrngui.hoc")

Default_Eleak = -65
membranecap = 0.64      	/* specific membrane capacitance in uF cm^-2 */
membraneresist = 120236 	/* specific membrane resistance in ohm cm^2 */
axialresist = 120	     	/* axial resistivity in ohm cm */

//Load the morphology file
xopen("PurkinjeMorphology_allbr.nrn")
forsec "axon" delete_section()	// Delete original axon and add AIS

objref g2, b2,c2, distrx, distry, cdistry, p

	forall {
		insert pas e_pas=Default_Eleak	/* Insert Leak everywhere */
	    insert hpkj	// Ih inserted everywhere
		insert ds
		insert pk
	}

    AIS {  g_pas=1/membraneresist Ra=axialresist cm=membranecap}
		forsec spinydend {g_pas=5.3/membraneresist Ra=axialresist cm=5.3*membranecap}
		forsec maindend {g_pas=1.2/membraneresist Ra=axialresist cm=1.2*membranecap}
		forsec "soma" { g_pas=1/membraneresist Ra=axialresist cm=membranecap}

		forsec maindend {insert cdp4N}
		forsec alldend {
				insert Kv3
				gbar_Kv3 = 0.1512
				vshift_Kv3 = 4
				insert newCaP
				pcabar_newCaP = 1.90e-4
				vshift_newCaP =-5
				insert CaT3_1
				pcabar_CaT3_1 = 2.7e-5
				insert mslo
				gbar_mslo = 0.21504
				insert SK2
				gkbar_SK2 = 3.6000e-04
				scal_SK2 = 1.0
				ghbar_hpkj = 0.00036
				insert Kv1
				gbar_Kv1 = 0.002//*2
				insert Kv4
				gbar_Kv4 = 0.0252
				insert Kv4s
				gbar_Kv4s = 0.015

		}
		forsec spinydend {
				insert cdp4Nsp
				gkbar_SK2 = 3.6000e-04
				scal_SK2 = 1.0
				gbar_Kv4 = 0.0264
				gbar_Kv4s = 0.015
				ghbar_hpkj = 3.2400e-04
				vshift_Kv4 = 0
				gbar_Kv1 = 0.001
				gbar_Kv3 = 0.1512
				vshift_Kv3 = 0
				pcabar_CaT3_1 = 1.0800e-04
				pcabar_newCaP = 7.6000e-04
				vshift_newCaP = -5
				scale_cdp4Nsp = 3.5
				gbar_mslo = 0.0896
				insert abBK
				gabkbar_abBK = 0.15
		}

// Include heterogeity for every branch of the PC:
		forsec branch1 {
		    insert cdp4Nsp
		    gkbar_SK2 = 3.6000e-04
		    scal_SK2 = 1.0
		    gbar_Kv4 = 0.0264
		    gbar_Kv4s = 0.015
		    ghbar_hpkj = 3.2400e-04
		    vshift_Kv4 = 0
		    gbar_Kv1 = 0.001
		    gbar_Kv3 = 0.1512
		    vshift_Kv3 = 0
		    pcabar_CaT3_1 = 1.0800e-04
		    pcabar_newCaP = 7.6000e-04
		    vshift_newCaP = -5
		    scale_cdp4Nsp = 3.5
		    gbar_mslo = 0.0896
		    insert abBK
		    gabkbar_abBK = 0.15
		}

		forsec branch2 {
		    insert cdp4Nsp
		    gkbar_SK2 = 3.6000e-04
		    scal_SK2 = 1.0
		    gbar_Kv4 = 0.0264
		    gbar_Kv4s = 0.015
		    ghbar_hpkj = 3.2400e-04
		    vshift_Kv4 = 0
		    gbar_Kv1 = 0.001
		    gbar_Kv3 = 0.1512
		    vshift_Kv3 = 0
		    pcabar_CaT3_1 = 1.0800e-04
		    pcabar_newCaP  = 7.6000e-04
		    vshift_newCaP = -5
		    scale_cdp4Nsp = 3.5
		    gbar_mslo = 0.0896
		    insert abBK
		    gabkbar_abBK = 0.15
		}

		forsec branch3 {
		    insert cdp4Nsp
		    gkbar_SK2 = 3.6000e-04
		    scal_SK2 = 1.0
		    gbar_Kv4 = 0.0264
		    gbar_Kv4s = 0.015
		    ghbar_hpkj = 3.2400e-04
		    vshift_Kv4 = 0
		    gbar_Kv1 = 0.001
		    gbar_Kv3 = 0.1512
		    vshift_Kv3 = 0
		    pcabar_CaT3_1 = 1.0800e-04
		    pcabar_newCaP = 7.6000e-04*1
		    vshift_newCaP = -5
		    scale_cdp4Nsp = 3.5
		    gbar_mslo = 0.0896
		    insert abBK
		    gabkbar_abBK = 0.15
		}

		forsec branch4 {
		    insert cdp4Nsp
		    gkbar_SK2 = 3.6000e-04
		    scal_SK2 = 1.0
		    gbar_Kv4 = 0.0264
		    gbar_Kv4s = 0.015
		    ghbar_hpkj = 3.2400e-04
		    vshift_Kv4 = 0
		    gbar_Kv1 = 0.001
		    gbar_Kv3 = 0.1512
		    vshift_Kv3 = 0
		    pcabar_CaT3_1 = 1.0800e-04
		    pcabar_newCaP = 7.6000e-04
		    vshift_newCaP = -5
		    scale_cdp4Nsp = 3.5
		    gbar_mslo = 0.0896
		    insert abBK
		    gabkbar_abBK = 0.15
		}

		forsec branch5 {
		    insert cdp4Nsp
		    gkbar_SK2 = 3.6000e-04
		    scal_SK2 = 1.0
		    gbar_Kv4 = 0.0264
		    gbar_Kv4s = 0.015
		    ghbar_hpkj = 3.2400e-04
		    vshift_Kv4 = 0
		    gbar_Kv1 = 0.001
		    gbar_Kv3 = 0.1512
		    vshift_Kv3 = 0
		    pcabar_CaT3_1 = 1.0800e-04
		    pcabar_newCaP = 7.6000e-04
		    vshift_newCaP = -5
		    scale_cdp4Nsp = 3.5
		    gbar_mslo = 0.0896
		    insert abBK
		    gabkbar_abBK = 0.15
		}

		forsec branch6 {
		    insert cdp4Nsp
		    gkbar_SK2 = 3.6000e-04
		    scal_SK2 = 1.0
		    gbar_Kv4 = 0.0264
		    gbar_Kv4s = 0.015
		    ghbar_hpkj = 3.2400e-04
		    vshift_Kv4 = 0
		    gbar_Kv1 = 0.001
		    gbar_Kv3 =  0.1512
		    vshift_Kv3 = 0
		    pcabar_CaT3_1 = 1.0800e-04
		    pcabar_newCaP = 7.6000e-04
		    vshift_newCaP = -5
		    scale_cdp4Nsp = 3.5
		    gbar_mslo = 0.0896
		    insert abBK
		    gabkbar_abBK = 0.15
		}

		forsec branch7 {
		    insert cdp4Nsp
		    gkbar_SK2 = 3.6000e-04
		    scal_SK2 = 1.0
		    gbar_Kv4 = 0.0264
		    gbar_Kv4s = 0.015
		    ghbar_hpkj = 3.2400e-04
		    vshift_Kv4 = 0
		    gbar_Kv1 = 0.001
		    gbar_Kv3 =  0.1512
		    vshift_Kv3 = 0
		    pcabar_CaT3_1 = 1.0800e-04
		    pcabar_newCaP = 7.6000e-04*2.8
		    vshift_newCaP = -5
		    scale_cdp4Nsp = 3.5
		    gbar_mslo = 0.0896
		    insert abBK
		    gabkbar_abBK = 0.15
		}

		forsec branch8 {
		    insert cdp4Nsp
		    gkbar_SK2 = 3.6000e-04
		    scal_SK2 = 1.0
		    gbar_Kv4 = 0.0264
		    gbar_Kv4s = 0.015
		    ghbar_hpkj = 3.2400e-04
		    vshift_Kv4 = 0
		    gbar_Kv1 = 0.001
		    gbar_Kv3 = 0.1512
		    vshift_Kv3 = 0
		    pcabar_CaT3_1 = 1.0800e-04
		    pcabar_newCaP  = 7.6000e-04*0.9
		    vshift_newCaP = -5
		    scale_cdp4Nsp = 3.5
		    gbar_mslo = 0.0896
		    insert abBK
		    gabkbar_abBK = 0.15
		}

		forsec branch9 {
		    insert cdp4Nsp
		    gkbar_SK2 = 3.6000e-04
		    scal_SK2 = 1.0
		    gbar_Kv4 = 0.0264
		    gbar_Kv4s = 0.015
		    ghbar_hpkj = 3.2400e-04
		    vshift_Kv4 = 0
		    gbar_Kv1 = 0.001
		    gbar_Kv3 = 0.1512
		    vshift_Kv3 = 0
		    pcabar_CaT3_1 = 1.0800e-04
		    pcabar_newCaP = 7.6000e-04
		    vshift_newCaP = -5
		    scale_cdp4Nsp = 3.5
		    gbar_mslo = 0.0896
		    insert abBK
		    gabkbar_abBK = 0.15
		}

		forsec branch10 {
		    insert cdp4Nsp
		    gkbar_SK2 = 3.6000e-04
		    scal_SK2 = 1.0
		    gbar_Kv4 = 0.0264
		    gbar_Kv4s = 0.015
		    ghbar_hpkj = 3.2400e-04
		    vshift_Kv4 = 0
		    gbar_Kv1 = 0.001
		    gbar_Kv3 = 0.1512
		    vshift_Kv3 = 0
		    pcabar_CaT3_1 = 1.0800e-04
		    pcabar_newCaP = 7.6000e-04
		    vshift_newCaP = -5
		    scale_cdp4Nsp = 3.5
		    gbar_mslo = 0.0896
		    insert abBK
		    gabkbar_abBK = 0.15
		}

		forsec branch11 {
		    insert cdp4Nsp
		    gkbar_SK2 = 3.6000e-04
		    scal_SK2 = 1.0
		    gbar_Kv4 = 0.0264
		    gbar_Kv4s = 0.015
		    ghbar_hpkj = 3.2400e-04
		    vshift_Kv4 = 0
		    gbar_Kv1 = 0.001
		    gbar_Kv3 =  0.1512
		    vshift_Kv3 = 0
		    pcabar_CaT3_1 = 1.0800e-04
		    pcabar_newCaP = 7.6000e-04
		    vshift_newCaP = -5
		    scale_cdp4Nsp = 3.5
		    gbar_mslo = 0.0896
		    insert abBK
		    gabkbar_abBK = 0.15
		}

		forsec branch12 {
		    insert cdp4Nsp
		    gkbar_SK2 = 3.6000e-04
		    scal_SK2 = 1.0
		    gbar_Kv4 = 0.0264
		    gbar_Kv4s = 0.015
		    ghbar_hpkj = 3.2400e-04
		    vshift_Kv4 = 0
		    gbar_Kv1 = 0.001
		    gbar_Kv3 = 0.1512
		    vshift_Kv3 = 0
		    pcabar_CaT3_1 = 1.0800e-04
		    pcabar_newCaP = 7.6000e-04*0.8
		    vshift_newCaP = -5
		    scale_cdp4Nsp = 3.5
		    gbar_mslo = 0.0896
		    insert abBK
		    gabkbar_abBK = 0.15
		}

		forsec branch13 {
		    insert cdp4Nsp
		    gkbar_SK2 = 3.6000e-04
		    scal_SK2 = 1.0
		    gbar_Kv4 = 0.0264
		    gbar_Kv4s = 0.015
		    ghbar_hpkj = 3.2400e-04
		    vshift_Kv4 = 0
		    gbar_Kv1 = 0.001
		    gbar_Kv3 = 0.1512
		    vshift_Kv3 = 0
		    pcabar_CaT3_1 = 1.0800e-04
		    pcabar_newCaP = 7.6000e-04
		    vshift_newCaP = -5
		    scale_cdp4Nsp = 3.5
		    gbar_mslo = 0.0896
		    insert abBK
		    gabkbar_abBK = 0.15
		}

		forsec branch14 {
		    insert cdp4Nsp
		    gkbar_SK2 = 3.6000e-04
		    scal_SK2 = 1.0
		    gbar_Kv4 = 0.0264
		    gbar_Kv4s = 0.015
		    ghbar_hpkj = 3.2400e-04
		    vshift_Kv4 = 0
		    gbar_Kv1 = 0.001
		    gbar_Kv3 =  0.1512
		    vshift_Kv3 = 0
		    pcabar_CaT3_1 = 1.0800e-04
		    pcabar_newCaP = 7.6000e-04
		    vshift_newCaP = -5
		    scale_cdp4Nsp = 3.5
		    gbar_mslo = 0.0896
		    insert abBK
		    gabkbar_abBK = 0.15
		}

		forsec branch15 {
		    insert cdp4Nsp
		    gkbar_SK2 = 3.6000e-04
		    scal_SK2 = 1.0
		    gbar_Kv4 = 0.0264
		    gbar_Kv4s = 0.015
		    ghbar_hpkj = 3.2400e-04
		    vshift_Kv4 = 0
		    gbar_Kv1 = 0.001
		    gbar_Kv3 =  0.1512
		    vshift_Kv3 = 0
		    pcabar_CaT3_1 = 1.0800e-04
		    pcabar_newCaP = 7.6000e-04
		    vshift_newCaP = -5
		    scale_cdp4Nsp = 3.5
		    gbar_mslo = 0.0896
		    insert abBK
		    gabkbar_abBK = 0.15
		}

		forsec branch16 {
		    insert cdp4Nsp
		    gkbar_SK2 = 3.6000e-04
		    scal_SK2 = 1.0
		    gbar_Kv4 = 0.0264
		    gbar_Kv4s = 0.015
		    ghbar_hpkj = 3.2400e-04
		    vshift_Kv4 = 0
		    gbar_Kv1 = 0.001
		    gbar_Kv3 = 0.1512
		    vshift_Kv3 = 0
		    pcabar_CaT3_1 = 1.0800e-04
		    pcabar_newCaP = 7.6000e-04
		    vshift_newCaP = -5
		    scale_cdp4Nsp = 3.5
		    gbar_mslo = 0.0896
		    insert abBK
		    gabkbar_abBK = 0.15
		}

		forsec branch17 {
		    insert cdp4Nsp
		    gkbar_SK2 = 3.6000e-04
		    scal_SK2 = 1.0
		    gbar_Kv4 = 0.0264
		    gbar_Kv4s = 0.015
		    ghbar_hpkj = 3.2400e-04
		    vshift_Kv4 = 0
		    gbar_Kv1 = 0.001
		    gbar_Kv3 = 0.1512
		    vshift_Kv3 = 0
		    pcabar_CaT3_1 = 1.0800e-04
		    pcabar_newCaP = 7.6000e-04
		    vshift_newCaP = -5
		    scale_cdp4Nsp = 3.5
		    gbar_mslo = 0.0896
		    insert abBK
		    gabkbar_abBK = 0.15
		}

		forsec branch18 {
		    insert cdp4Nsp
		    gkbar_SK2 = 3.6000e-04
		    scal_SK2 = 1.0
		    gbar_Kv4 = 0.0264*1.0
		    gbar_Kv4s = 0.015
		    ghbar_hpkj = 3.2400e-04
		    vshift_Kv4 = 0
		    gbar_Kv1 = 0.001
		    gbar_Kv3 = 0.1512
		    vshift_Kv3 = 0
		    pcabar_CaT3_1 = 1.0800e-04
		    pcabar_newCaP = 7.6000e-04
		    vshift_newCaP = -5
		    scale_cdp4Nsp = 3.5
		    gbar_mslo = 0.0896
		    insert abBK
		    gabkbar_abBK = 0.15
		}

		forsec branch19 {
		    insert cdp4Nsp
		    gkbar_SK2 = 3.6000e-04
		    scal_SK2 = 1.0
		    gbar_Kv4 = 0.0264*1.0
		    gbar_Kv4s  = 0.015
		    ghbar_hpkj = 3.2400e-04
		    vshift_Kv4 = 0
		    gbar_Kv1 = 0.001
		    gbar_Kv3 = 0.1512
		    vshift_Kv3 = 0
		    pcabar_CaT3_1 = 1.0800e-04
		    pcabar_newCaP = 7.6000e-04
		    vshift_newCaP = -5
		    scale_cdp4Nsp = 3.5
		    gbar_mslo = 0.0896
		    insert abBK
		    gabkbar_abBK = 0.15
		}

		forsec branch20 {
		    insert cdp4Nsp
		    gkbar_SK2 = 3.6000e-04
		    scal_SK2 = 1.0
		    gbar_Kv4 = 0.0264*1.0
		    gbar_Kv4s = 0.015
		    ghbar_hpkj = 3.2400e-04
		    vshift_Kv4 = 0
		    gbar_Kv1 = 0.001
		    gbar_Kv3 =  0.1512
		    vshift_Kv3 = 0
		    pcabar_CaT3_1 = 1.0800e-04
		    pcabar_newCaP = 7.6000e-04
		    vshift_newCaP = -5
		    scale_cdp4Nsp = 3.5
		    gbar_mslo = 0.0896
		    insert abBK
		    gabkbar_abBK = 0.15
		}

		forsec branch21 {
		    insert cdp4Nsp
		    gkbar_SK2 = 3.6000e-04
		    scal_SK2 = 1.0
		    gbar_Kv4 = 0.0264
		    gbar_Kv4s = 0.015
		    ghbar_hpkj = 3.2400e-04
		    vshift_Kv4 = 0
		    gbar_Kv1 = 0.001
		    gbar_Kv3 =  0.1512
		    vshift_Kv3 = 0
		    pcabar_CaT3_1 = 1.0800e-04
		    pcabar_newCaP  =7.6000e-04*0.8
		    vshift_newCaP = -5
		    scale_cdp4Nsp = 3.5
		    gbar_mslo = 0.0896
		    insert abBK
		    gabkbar_abBK = 0.15
		}

		forsec branch22 {
		    insert cdp4Nsp
		    gkbar_SK2 = 3.6000e-04
		    scal_SK2 = 1.0
		    gbar_Kv4 = 0.0264
		    gbar_Kv4s = 0.015
		    ghbar_hpkj = 3.2400e-04
		    vshift_Kv4 = 0
		    gbar_Kv1 = 0.001
		    gbar_Kv3 = 0.1512
		    vshift_Kv3 = 0
		    pcabar_CaT3_1 = 1.0800e-04
		    pcabar_newCaP = 7.6000e-04*0.8
		    vshift_newCaP = -5
		    scale_cdp4Nsp = 3.5
		    gbar_mslo = 0.0896
		    insert abBK
		    gabkbar_abBK = 0.15
		}

access somaA
somaA distance(0,0.5)

forsec "soma" {
	 // INa & IK's at soma only
	//for soma, I make it fire under dissociated condition. The firing properties is comparable with Bean's data
			insert naRsg		// INa + resurgent
			gbar_naRsg = 0.03168
			vshifta_naRsg = 0
			vshiftk_naRsg = 0
			vshifti_naRsg = -5

			insert nap
			gbar_nap =  1.4000e-04
			insert pk
			ena = 63
			ghbar_hpkj = 1.0800e-04 // Default = 1

			insert cdp20N_FD2
			insert Kv3
			gbar_Kv3 = 1.8
			vshift_Kv3 = 4

			insert Kv1
			gbar_Kv1 = 0
			insert newCaP
			pcabar_newCaP = 1.9000e-04
			kt_newCaP = 1
			vshift_newCaP = -5
			insert mslo
			gbar_mslo = 0.8736
			insert abBK
			gabkbar_abBK = 0.3
			insert SK2
			gkbar_SK2 = 0.0075
	}

	AIS {
				insert naRsg		// INa + resurgent
				gbar_naRsg = 0.56
				vshifta_naRsg = 15
				vshiftk_naRsg = 5
				vshifti_naRsg = -5

				insert nap
				gbar_nap = 0.0023

				insert CaT3_1
				pcabar_CaT3_1 =  1.2800e-04
				ena = 63
				ghbar_hpkj =  1.0800e-04


				insert cdpAIS
				insert Kv3
				gbar_Kv3 = 115.2
				vshift_Kv3 = 4
				insert Kv1
				gbar_Kv1 = 0

				insert newCaP
				pcabar_newCaP = 0.00228
				kt_newCaP = 1
				vshift_newCaP = -5

				insert mslo
				gbar_mslo = 6
				insert abBK
				gabkbar_abBK = 1.05
				insert SK2
				gkbar_SK2 = 0.0278
		}


objref patch_site
patch_site = new List()

ip = 0
forsec patch_list {
patch_site.append(new SectionRef())
}

celsius = 34
dt = 0.02
steps_per_ms = 1/dt
dtsim = 0.02

objref g2, b2,c2, distrm, distrd

xopen ("electrode.hoc")
xopen("distri.hoc")	//voltage spatial distribution

proc clamp_cc() {

	somaA {
		stim1.del = 0
		stim1.dur = 1000000
		stim1.amp = $1
	}
}

v_init = -70

// Define the length of our simulation:
tstop = 600

objref scalefile
scalefile=new File()

xopen("distri_synapse.hoc")

objref pc
pc = new ParallelContext()

//function farmed out to slave nodes
func distscale() {local key, errval, cu_id, fr_id, dend_id localobj parvec, returnlist
	key = $1
	cu_id = int($1/100000000)
	fr_id = int(($1 - cu_id*100000000)/1000000)
	site_id = int(($1 - cu_id*100000000-fr_id*1000000)/10000)
	dend_id = int(($1 - cu_id*100000000-fr_id*1000000-site_id*10000)/100)
	trial_id = $1 - cu_id*100000000-fr_id*1000000-site_id*10000-dend_id*100
	returnlist = new List()
	returnlist = calc_EPSP_single(cu_id,fr_id,site_id,dend_id,trial_id)

	pc.pack(returnlist.o(0))
	pc.pack(returnlist.o(1))
	pc.pack(returnlist.o(2))
	pc.pack(returnlist.o(3))
	pc.pack(returnlist.o(4))
	pc.pack(returnlist.o(5))
	pc.pack(returnlist.o(6))
	pc.pack(returnlist.o(7))
	pc.pack(returnlist.o(8))
	pc.pack(returnlist.o(9))
	pc.pack(returnlist.o(10))
	pc.pack(returnlist.o(11))
	pc.pack(returnlist.o(12))
	pc.pack(returnlist.o(13))
	pc.pack(returnlist.o(14))
	pc.pack(returnlist.o(15))
	pc.pack(returnlist.o(16))
	pc.pack(returnlist.o(17))
	pc.pack(returnlist.o(18))
	pc.pack(returnlist.o(19))
	pc.post(key)
	return key
}

objref aSynapseList[11]
obfunc calc_EPSP_single() {localobj outlist, currecord, integ_soma, br1,br1_dist,br1_prox,br2_dist,br2_prox,br3_dist,br3_prox,br4_dist,br4_prox,br5_dist,br5_dist2,br5_prox,br5_prox2,br6_dist,br6_prox,br7_dist,br7_prox,br8_dist,br8_prox,br8_mid
	//function to calculate the max deflection due to a single synapse
	cu_id = $1
	fr_id = $2
	siteval = $3
	synval= $4
	tr_id = $5
	curr = -0.4+(cu_id-1)*0.2

	nlist = fr_id	// in fact nlist can be multiple, make synapses firing at bursting

	Npf = 2+(synval-1)*2 // Number of PF

for i = 1,nlist {aSynapseList[i-1] = new List() }	// every time this will be initialized.
randomseed = cu_id*100000000+fr_id*1000000+siteval*10000 +synval*100 + tr_id

br1 = new SectionList()
patch_site.o(siteval-1).sec br1.subtree()
br1.remove(cf)
br1.unique()

if (siteval==12) {
	br1 = new SectionList()
	forsec "dendA1_001101100*" br1.append()
}
if (i==13) {
	br1 = new SectionList()
	forsec "dendA1_0011011100*" br1.append()
	forsec "dendA1_0011011101*" br1.append()
}

for i = 1,nlist {aSynapseList[i-1] = distSyns(Npf,br1,randomseed)}

for i = 1,nlist {
	for ii=1,aSynapseList[i-1].count() {
  	  aSynapseList[i-1].object(ii-1).onset = 399 //define start time of the PF activation
  	  aSynapseList[i-1].object(ii-1).tau0 = 0.3
  	  aSynapseList[i-1].object(ii-1).tau1 = 3
   	 aSynapseList[i-1].object(ii-1).e = 0
   	 aSynapseList[i-1].object(ii-1).gmax = 0.5e-3//
	}
}

// Record vectors containing results for the soma and for different distal and proximal points on the branches situated on the left main division of the tree:
outlist=new List()
integ_soma = new Vector()
integ_soma.record(&somaA.v(0.5))
br1_dist = new Vector()
br1_dist.record(&dendA1_00001101.v(0.05))
br1_prox = new Vector()
br1_prox.record(&dendA1_000.v(0.05))
br2_dist = new Vector()
br2_dist.record(&dendA1_001000111.v(0.05))
br2_prox = new Vector()
br2_prox.record(&dendA1_00100.v(0.05))
br3_dist = new Vector()
br3_dist.record(&dendA1_0010101100.v(0.05))
br3_prox = new Vector()
br3_prox.record(&dendA1_001010.v(0.05))
br4_dist = new Vector()
br4_dist.record(&dendA1_00101101001111.v(0.05))
br4_prox = new Vector()
br4_prox.record(&dendA1_00101100.v(0.05))
br5_dist = new Vector()
br5_dist.record(&dendA1_00101110101000110.v(0.05))
br5_dist2 = new Vector()
br5_dist2.record(&dendA1_0010111011000.v(0.05))
br5_prox = new Vector()
br5_prox.record(&dendA1_001011101111.v(0.05))
br5_prox2 = new Vector()
br5_prox2.record(&dendA1_0010111010.v(0.05))
br6_dist = new Vector()
br6_dist.record(&dendA1_0010111100111.v(0.05))
br6_prox = new Vector()
br6_prox.record(&dendA1_0010111100.v(0.05))
br7_dist = new Vector()
br7_dist.record(&dendA1_0010111101011001111.v(0.05))
br7_prox = new Vector()
br7_prox.record(&dendA1_00101111010.v(0.05))
br8_dist = new Vector()
br8_dist.record(&dendA1_001011110110010101.v(0.05))
br8_prox = new Vector()
br8_prox.record(&dendA1_001011110110.v(0.05))
br8_mid = new Vector()
br8_mid.record(&dendA1_00101111011001.v(0.05))

clamp_cc(curr)
finitialize(v_init)
continuerun(tstop)

outlist.append(integ_soma)
outlist.append(br1_dist)
outlist.append(br1_prox)
outlist.append(br2_dist)
outlist.append(br2_prox)
outlist.append(br3_dist)
outlist.append(br3_prox)
outlist.append(br4_dist)
outlist.append(br4_prox)
outlist.append(br5_dist)
outlist.append(br5_dist2)
outlist.append(br5_prox)
outlist.append(br5_prox2)
outlist.append(br6_dist)
outlist.append(br6_prox)
outlist.append(br7_dist)
outlist.append(br7_prox)
outlist.append(br8_dist)
outlist.append(br8_prox)
outlist.append(br8_mid)


return outlist
}

pc.runworker()

objref somavec,br1_distvec, br1_proxvec, br2_distvec, br2_proxvec,br3_distvec, br3_proxvec,br4_distvec, br4_proxvec, br4_distvec, br4_proxvec,br5_distvec, br5_distvec2,br5_proxvec,br5_proxvec2,br6_distvec,br6_proxvec,br7_distvec,br7_proxvec,br8_distvec,br8_proxvec,br8_midvec

somavec = new Vector()
br1_distvec = new Vector()
br1_proxvec = new Vector()
br2_distvec = new Vector()
br2_proxvec = new Vector()
br3_distvec = new Vector()
br3_proxvec = new Vector()
br4_distvec = new Vector()
br4_proxvec = new Vector()
br5_distvec = new Vector()
br5_distvec2 = new Vector()
br5_proxvec = new Vector()
br5_proxvec2 = new Vector()
br6_distvec = new Vector()
br6_proxvec = new Vector()
br7_distvec = new Vector()
br7_proxvec = new Vector()
br8_distvec = new Vector()
br8_proxvec = new Vector()
br8_midvec = new Vector()



strdef tmpstr
strdef outDir
strdef cmd
objref outfile
outfile = new File()

proc calcEPSPs() {
	sprint(outDir,"Results")
	sprint(cmd, "system(\"mkdir -p %s\")",outDir)
	execute(cmd)
	somaA distance(0,0.5)
	for cu = 3, $1 {
		for freq = 1, $2 {
			for site = $3, $3 {
				for m = $4, $4 {
					for nt =$5 ,$5 {
						mmtag=cu*100000000 + freq*1000000 + site*10000+m*100 + nt
						pc.submit("distscale",mmtag)	//send out the error calculations
					}
				}
			}
		}
	}
	//collect error values
	while (pc.working()) {
		key = pc.retval()	//retrieve the tag
		pc.look_take(key)	//remove the tag/job from the bulletin
		somavec = pc.upkvec()	//unpack the error value associated with the tag
		br1_distvec = pc.upkvec()
		br1_proxvec = pc.upkvec()
		br2_distvec = pc.upkvec()
		br2_proxvec = pc.upkvec()
		br3_distvec = pc.upkvec()
		br3_proxvec = pc.upkvec()
		br4_distvec = pc.upkvec()
		br4_proxvec = pc.upkvec()
		br5_distvec = pc.upkvec()
		br5_distvec2 = pc.upkvec()
		br5_proxvec = pc.upkvec()
		br5_proxvec2 = pc.upkvec()
		br6_distvec = pc.upkvec()
		br6_proxvec = pc.upkvec()
		br7_distvec = pc.upkvec()
		br7_proxvec = pc.upkvec()
		br8_distvec = pc.upkvec()
		br8_proxvec = pc.upkvec()
		br8_midvec = pc.upkvec()

		print "received key ",key
		cuno = int(key/100000000)
		frno = int((key- cuno*100000000)/1000000)
		siteno = int((key - cuno*100000000 - frno*1000000)/10000)
		synno= int((key - cuno*100000000 - frno*1000000 - siteno*10000)/100)
		trno  = key - cuno*100000000-frno*1000000-siteno*10000-synno*100

		sprint(tmpstr,"%s/%03d_%03d_%03d_%03d_%03d_vsoma.dat",outDir,cuno,frno,siteno,synno,trno)
		outfile.wopen(tmpstr)
		somavec.printf(outfile)
		outfile.close()

		sprint(tmpstr,"%s/%03d_%03d_%03d_%03d_%03d_vbr1_dist.dat",outDir,cuno,frno,siteno,synno,trno)
		outfile.wopen(tmpstr)
		br1_distvec.printf(outfile)
		outfile.close()

		sprint(tmpstr,"%s/%03d_%03d_%03d_%03d_%03d_vbr1_prox.dat",outDir,cuno,frno,siteno,synno,trno)
		outfile.wopen(tmpstr)
		br1_proxvec.printf(outfile)
		outfile.close()

		sprint(tmpstr,"%s/%03d_%03d_%03d_%03d_%03d_vbr2_dist.dat",outDir,cuno,frno,siteno,synno,trno)
		outfile.wopen(tmpstr)
		br2_distvec.printf(outfile)
		outfile.close()

		sprint(tmpstr,"%s/%03d_%03d_%03d_%03d_%03d_vbr2_prox.dat",outDir,cuno,frno,siteno,synno,trno)
		outfile.wopen(tmpstr)
		br2_proxvec.printf(outfile)
		outfile.close()

		sprint(tmpstr,"%s/%03d_%03d_%03d_%03d_%03d_vbr3_dist.dat",outDir,cuno,frno,siteno,synno,trno)
		outfile.wopen(tmpstr)
		br3_distvec.printf(outfile)
		outfile.close()

		sprint(tmpstr,"%s/%03d_%03d_%03d_%03d_%03d_vbr3_prox.dat",outDir,cuno,frno,siteno,synno,trno)
		outfile.wopen(tmpstr)
		br3_proxvec.printf(outfile)
		outfile.close()

		sprint(tmpstr,"%s/%03d_%03d_%03d_%03d_%03d_vbr4_dist.dat",outDir,cuno,frno,siteno,synno,trno)
		outfile.wopen(tmpstr)
		br4_distvec.printf(outfile)
		outfile.close()

		sprint(tmpstr,"%s/%03d_%03d_%03d_%03d_%03d_vbr4_prox.dat",outDir,cuno,frno,siteno,synno,trno)
		outfile.wopen(tmpstr)
		br4_proxvec.printf(outfile)
		outfile.close()

		sprint(tmpstr,"%s/%03d_%03d_%03d_%03d_%03d_vbr5_dist.dat",outDir,cuno,frno,siteno,synno,trno)
		outfile.wopen(tmpstr)
		br5_distvec.printf(outfile)
		outfile.close()

		sprint(tmpstr,"%s/%03d_%03d_%03d_%03d_%03d_vbr5_dist2.dat",outDir,cuno,frno,siteno,synno,trno)
		outfile.wopen(tmpstr)
		br5_distvec2.printf(outfile)
		outfile.close()

		sprint(tmpstr,"%s/%03d_%03d_%03d_%03d_%03d_vbr5_prox.dat",outDir,cuno,frno,siteno,synno,trno)
		outfile.wopen(tmpstr)
		br5_proxvec.printf(outfile)
		outfile.close()

		sprint(tmpstr,"%s/%03d_%03d_%03d_%03d_%03d_vbr5_prox2.dat",outDir,cuno,frno,siteno,synno,trno)
		outfile.wopen(tmpstr)
		br5_proxvec2.printf(outfile)
		outfile.close()

		sprint(tmpstr,"%s/%03d_%03d_%03d_%03d_%03d_vbr6_dist.dat",outDir,cuno,frno,siteno,synno,trno)
		outfile.wopen(tmpstr)
		br6_distvec.printf(outfile)
		outfile.close()

		sprint(tmpstr,"%s/%03d_%03d_%03d_%03d_%03d_vbr6_prox.dat",outDir,cuno,frno,siteno,synno,trno)
		outfile.wopen(tmpstr)
		br6_proxvec.printf(outfile)
		outfile.close()

		sprint(tmpstr,"%s/%03d_%03d_%03d_%03d_%03d_vbr7_dist.dat",outDir,cuno,frno,siteno,synno,trno)
		outfile.wopen(tmpstr)
		br7_distvec.printf(outfile)
		outfile.close()

		sprint(tmpstr,"%s/%03d_%03d_%03d_%03d_%03d_vbr7_prox.dat",outDir,cuno,frno,siteno,synno,trno)
		outfile.wopen(tmpstr)
		br7_proxvec.printf(outfile)
		outfile.close()

		sprint(tmpstr,"%s/%03d_%03d_%03d_%03d_%03d_vbr8_dist.dat",outDir,cuno,frno,siteno,synno,trno)
		outfile.wopen(tmpstr)
		br8_distvec.printf(outfile)
		outfile.close()

		sprint(tmpstr,"%s/%03d_%03d_%03d_%03d_%03d_vbr8_prox.dat",outDir,cuno,frno,siteno,synno,trno)
		outfile.wopen(tmpstr)
		br8_proxvec.printf(outfile)
		outfile.close()

		sprint(tmpstr,"%s/%03d_%03d_%03d_%03d_%03d_vbr8_mid.dat",outDir,cuno,frno,siteno,synno,trno)
		outfile.wopen(tmpstr)
		br8_midvec.printf(outfile)
		outfile.close()
	}
}

calcEPSPs(3,1,7,21,1)
// This codes reproduces Figure 4G-I from Cirtala2024
// Therefore only the case corresponding to br7 - for 42PF (branch 7 PF threshold value)
// If you wish to run it for the entire range of PF please change the following line inside of the proc calcESP:
// for m = $4, $4  ---> for m = 1, $4
// and instead of 21 in calcEPSPs, please put the maximum number of PF you want (e.g. 75 corresponds to 150PF)