////////////////////
//Model Parameters//
////////////////////

////////////////////
//MF->GrC synapses//
////////////////////
distances=5.5556 //distances between GrCs
CONNECTION_WEIGHT_NMDA_MFGrC=0.2244e-3 //synaptic weight
CONNECTION_WEIGHT_AMPA_MFGrC=2.6e-3
CONNECTION_PROB_MFGrC=1  
CONNECTION_THRESHOLD_MFGrC=-20 
CONNECTION_DELAY_MFGrC=2e-3  
diameter_MFGrC=3.5*distances      
k1_MFGrC=0.01  //decay parameter

////////////////////
//MF->GoC synapses//
////////////////////
CONNECTION_WEIGHT_AMPA_MFGoC=2.6e-3  
CONNECTION_PROB_MFGoC=0.3
CONNECTION_THRESHOLD_MFGoC=-20 
CONNECTION_DELAY_MFGoC=2e-3   
diameter_MFGoC=52.5*distances
k1_MFGoC=0.01  

/////////////////////
//GrC->GoC synapses//
//(PF synapses)    //
/////////////////////
CONNECTION_WEIGHT_AMPA_GrCGoC=2.6e-3  
CONNECTION_PROB_GrCGoC=0.5
CONNECTION_THRESHOLD_GrCGoC_=-20 
CONNECTION_DELAY_GrCGoC=2e-3   
CONNECTION_DELAY_AF=300e-3  
k1_GrCGoC=0  

////////////////////
//GoC>GrC synapses//
////////////////////
CONNECTION_WEIGHT_GABAa1_GoCGrC=1.4e-3  
CONNECTION_WEIGHT_GABAa2_GoCGrC=0*1.4e-3  
CONNECTION_PROB_GoCGrC=1
CONNECTION_THRESHOLD_GoCGrC=-20 
CONNECTION_DELAY_GoCGrC=2e-3   
diameter_GoCGrC=10.5*distances 
k1_GoCGrC=0.01

/////////////////////////
//GoC>GoC Gap junctions//
////////////////////////

Ggap=1.5e-3  
CONNECTION_PROB_GoCGoC=1
k1_GoCGoC =0.03

//////////////////////
//Network Layers    //
//////////////////////

//load cells:
load_file("granule_template.hoc") //GrC
load_file("MF_template.hoc") //MF
load_file("Golgi_template.hoc") //GoC

//////////////////////////
//MF (Mossy Fibre) Layer//
//////////////////////////

scale=3  
objref prob        
prob=new Random()
NETDIM_Xa=90
NETDIM_Ya=10
NETDIM_Za=1
Dxa=scale*distances  
Dya=scale*distances
Dza=1
xa=0
ya=0
za=0

NumberCellsa=NETDIM_Xa*NETDIM_Ya*NETDIM_Za

variabilitya=20 
percentagea=variabilitya/100  
min=-Dxa*percentagea 
max=Dxa*percentagea

objref fiber[NumberCellsa]	//define objects for cells
l=0
for i=0, NETDIM_Xa-1 {
 for j=0, NETDIM_Ya-1 {
  for k=0, NETDIM_Za-1 {
  xprob=prob.uniform(min,max)
  yprob=prob.uniform(min,max)
 
   fiber[l]=new fibre(xprob+xa+i*Dxa,yprob+ya+j*Dya,za+k*Dza)  //create objects for cells
   fiber[l].Presynapse fiber[l].nclist.append(new NetCon(fiber[l].StimTrigger,fiber[l].syn,0.5,0,0.01)) 

   l=l+1
  }
 }
}

/////////////////////////////////////
//GrC (Granule Cell) Layer         //
/////////////////////////////////////

NETDIM_Xb= scale*NETDIM_Xa
NETDIM_Yb= scale*NETDIM_Ya
NETDIM_Zb=1
Dxb=distances	
Dyb=distances
Dzb=1
xb=0				
yb=0
zb=150

NumberCellsb=NETDIM_Xb*NETDIM_Yb*NETDIM_Zb

variabilityb=20 
percentageb=variabilityb/100  
min=-Dxb*percentageb 
max=Dxb*percentageb


objref GrC[NumberCellsb]	//define objects for cells
l=0
for i=0, NETDIM_Xb-1 {
 for j=0, NETDIM_Yb-1 {
  for k=0, NETDIM_Zb-1 {
  xprob=prob.uniform(min,max)
  yprob=prob.uniform(min,max)
  zprob=prob.uniform(min,0)
  GrC[l]=new grc(xprob+xb+i*Dxb,yprob+yb+j*Dyb,zprob+zb+k*Dzb)  //create objects for cells

 //soma diameter randomization
 diameter=GrC[l].soma.diam
 bprob=prob.uniform(-diameter*(variabilityb/100), diameter*(variabilityb/100))
  GrC[l].soma.diam= diameter +bprob

 //v_leak randomization
 v_leak=-16.5 
 dprob=prob.uniform(-v_leak*(variabilityb/100),v_leak*(variabilityb/100))
 GrC[l].soma.el_GRC_LKG1=v_leak+dprob
 l=l+1
  }
 }
}

//////////////////////////
//GoC (Golgi Cell) Layer//
//////////////////////////

scaled=2  
NETDIM_Xd=NETDIM_Xa/scaled
NETDIM_Yd=NETDIM_Ya/scaled
NETDIM_Zd=1

Dxd=scale*distances*scaled	
Dyd=scale*distances*scaled
Dzd=1
xd=0				
yd=0
zd=200

NumberCellsd=NETDIM_Xd*NETDIM_Yd*NETDIM_Zd

variabilityd=20  
percentaged= variabilityd/100  
min=-Dxd* percentaged 
max=Dxd* percentaged

Dxd=Dxd-2*(max/(NETDIM_Xd-1))
Dyd=Dyd-2*(max/(NETDIM_Xd-1))

xd=xd+max
yd=yd+max

objref GoC[NumberCellsd]	//define objects for cells
l=0
for i=0, NETDIM_Xd-1 {
 for j=0, NETDIM_Yd-1 {
  for k=0, NETDIM_Zd-1 {
  xprob=prob.uniform(min,max)
  yprob=prob.uniform(min,max)
  GoC[l]=new Goc(xprob+xd+i*Dxd,yprob+yd+j*Dyd,zd+k*Dzd)  //create objects for cells
  
v_leak=-55 
dprob=prob.uniform(-v_leak*(variabilityd/100),v_leak*(variabilityd/100))

GoC[l].soma.el_Golgi_lkg=v_leak+dprob
GoC[l].dend[0].el_Golgi_lkg=v_leak+dprob
GoC[l].dend[1].el_Golgi_lkg=v_leak+dprob
GoC[l].dend[2].el_Golgi_lkg=v_leak+dprob
//GoC[l].axon.el_Golgi_lkg=v_leak+dprob
GoC[l].soma.ko(0.5)=5
GoC[l].soma.ki(0.5)=140
GoC[l].soma.nao(0.5)=145
GoC[l].soma.nai(0.5)=5

  diameter =GoC[l].soma.diam
  dprob=prob.uniform(-diameter*(variabilityd/100), diameter*(variabilityd/100))
  GoC[l].soma.diam= diameter+dprob
   l=l+1
  }
 }
}

////////////////////////
//CONNECT CELL OBJECTS//
////////////////////////

objref prob
prob=new Random()

////////////////////
//MF->GrC synapses//
////////////////////

for i=0, NumberCellsa-1{
 for j=0, NumberCellsb-1 {
  aprob=prob.uniform(0,1)
  if ((aprob<= CONNECTION_PROB_MFGrC)) { 
    access fiber[i].Presynapse
    xi=x3d(0)
    yi=y3d(0)
    zi=0   
 
    access GrC[j].soma
    xj=x3d(0)
    yj=y3d(0)
    zj=0   
 
    adist =((xi-xj)^2+(yi-yj)^2+(zi-zj)^2)^(1/2)             			
    if (adist< diameter_MFGrC) {
   fiber[i].Presynapse GrC[j].nclist.append(new NetCon(&fiber[i].Presynapse.v(0.5),GrC[j].synAMPA, CONNECTION_THRESHOLD_MFGrC, CONNECTION_DELAY_MFGrC*adist,CONNECTION_WEIGHT_AMPA_MFGrC*(E^(-k1_MFGrC*adist))))
   fiber[i].Presynapse GrC[j].nclist.append(new NetCon(&fiber[i].Presynapse.v(0.5),GrC[j].synNMDA, CONNECTION_THRESHOLD_MFGrC, CONNECTION_DELAY_MFGrC*adist,CONNECTION_WEIGHT_NMDA_MFGrC*(E^(-k1_MFGrC*adist))))  
     }
  }
 }
}

////////////////////
//MF->GoC synapses//
////////////////////

for i=0, NumberCellsa-1{
 for j=0, NumberCellsd-1 {
  aprob=prob.uniform(0,1)
  if ((aprob<= CONNECTION_PROB_MFGoC)) { 
    access fiber[i].Presynapse
    xi=x3d(0)
    yi=y3d(0)
    zi=0  
 
    access  GoC[j].soma
    xj=x3d(0)
    yj=y3d(0)
    zj=0   
 
    adist =((xi-xj)^2+(yi-yj)^2+(zi-zj)^2)^(1/2)                 			
    if (adist< diameter_MFGoC) {  
  fiber[i].Presynapse GoC[j].nclist.append(new NetCon(&fiber[i].Presynapse.v(0.5),GoC[j].synAMPA1, CONNECTION_THRESHOLD_MFGoC, CONNECTION_DELAY_MFGoC*adist,CONNECTION_WEIGHT_AMPA_MFGoC*(E^(-k1_MFGoC*adist))))
     }
  }
 }
}

/////////////////////
//GrC->GoC synapses//
//(PF synapses)    //
/////////////////////

for i=0, NumberCellsb-1{
 for j=0, NumberCellsd-1 {
  aprob=prob.uniform(0,1)
  //aprob=1
  if ((aprob<= CONNECTION_PROB_GrCGoC)) {  
    access GrC[i].soma 
    xi=x3d(0)
    yi=y3d(0)
    zi=0   
    access GoC[j].soma 
    xj=x3d(0)
    yj=y3d(0)
    zj=0  
  
   adist =((xi-xj)^2+(yi-yj)^2+(zi-zj)^2)^(1/2)               			
   if (abs(xi-xj)<=2500 && abs(yi-yj)<=(2.6+100)) {
   GrC[i].soma GoC[j].nclist.append(new NetCon(&GrC[i].soma.v(0.5),GoC[j].synAMPA2, CONNECTION_THRESHOLD_GrCGoC_, CONNECTION_DELAY_GrCGoC*adist+CONNECTION_DELAY_AF,CONNECTION_WEIGHT_AMPA_GrCGoC*(E^(-k1_GrCGoC*adist))))
     }
  }
 }
}

////////////////////
//GoC>GrC synapses//
////////////////////

for i=0, NumberCellsd-1{
 for j=0, NumberCellsb-1 {
  aprob=prob.uniform(0,1)
  //aprob=1
  if ((aprob<= CONNECTION_PROB_GoCGrC)) { 
    access  GoC[i].soma  //GoC
    xi=x3d(0)
    yi=y3d(0)
    zi=0   
 
    access GrC[j].soma 
    xj=x3d(0)
    yj=y3d(0)
    zj=0   
 
    adist =((xi-xj)^2+(yi-yj)^2+(zi-zj)^2)^(1/2)             			
   if (adist< diameter_GoCGrC) {  
   GoC[i].soma GrC[j].nclist.append(new NetCon(&GoC[i].soma.v(0.5),GrC[j].synGABAa1, CONNECTION_THRESHOLD_GoCGrC, CONNECTION_DELAY_GoCGrC*adist,CONNECTION_WEIGHT_GABAa1_GoCGrC*(E^(-k1_GoCGrC*adist))))
   GoC[i].soma GrC[j].nclist.append(new NetCon(&GoC[i].soma.v(0.5),GrC[j].synGABAa2, CONNECTION_THRESHOLD_GoCGrC, CONNECTION_DELAY_GoCGrC*adist,CONNECTION_WEIGHT_GABAa2_GoCGrC*(E^(-k1_GoCGrC*adist))))
     }
  }
 }
}

/////////////////////////
//GoC>GoC Gap junctions//
////////////////////////

objref gpre[255*255],gpost[255*255]
objref index
index=new Matrix(255,255)

proc gapjunctions() {
 
contador=0
for i=0,  NumberCellsd-1{
for j=0, NumberCellsd-1 {
index.x[i][j]=contador
contador=contador+1
}
}

variabilitye=60  
percentagee= variabilitye/100 
min=-$3* percentagee 
max=$3* percentagee

for i=0, $4-1{
for j=i+1, $4-1 {
  aprob=1
  if ((aprob<= $1) && (i!=j)) {  
    access GoC[i].soma 
    xi=x3d(0)
    yi=y3d(0)
    zi=0   
 
    access GoC[j].soma 
    xj=x3d(0)
    yj=y3d(0)
    zj=0  
 
    adist =((xi-xj)^2+(yi-yj)^2+(zi-zj)^2)^(1/2)          			

GoC[i].soma gpre[index.getval(i,j)]=new Gap(0.5)
GoC[j].soma gpost[index.getval(i,j)]=new Gap(0.5)
setpointer gpre[index.getval(i,j)].vnb, GoC[j].soma.v(0.5)
setpointer gpost[index.getval(i,j)].vnb, GoC[i].soma.v(0.5)

gprob=prob.uniform(min,max) 

gpre[index.getval(i,j)].g=($3+gprob)*E^(-$2*adist)
gpost[index.getval(i,j)].g=($3+gprob)*E^(-$2*adist)
  }
 }
}
}

if (GJ==1) {
 gapjunctions(CONNECTION_PROB_GoCGoC, k1_GoCGoC, Ggap, NumberCellsd)
}


/////////////////////////////
//MF stimulation procedures//
////////////////////////////

//clear stim buffers
proc cleanMF() {     
for(i=0;i<=NumberCellsa-1;i=i+1) {
   fiber[i].StimTrigger.start=0
   fiber[i].StimTrigger.interval=0
   fiber[i].StimTrigger.number=0
   fiber[i].StimTrigger.noise=0
}
}

proc stimulus() {

NETIN_SP_START=$1 
MFfreq=$2
NETIN_SP_NUMBER=50*$3 
SEED=$4
NETIN_SP_INTERVAL=1000/MFfreq 
NETIN_SP_NOISE=1

for i = 0, NumberCellsa-1 {
   fiber[i].StimTrigger.start=NETIN_SP_START
   fiber[i].StimTrigger.interval=NETIN_SP_INTERVAL
   fiber[i].StimTrigger.number=NETIN_SP_NUMBER
   fiber[i].StimTrigger.noise=NETIN_SP_NOISE
   fiber[i].StimTrigger.seed(SEED)
}
}

////////////////////////////////////////
//Setting up for Saving Spiking Events//
///////////////////////////////////////

//generate matrix and send the data to them
objref GoCMatrixS, GrCMatrixS, MFMatrixS 
GoCMatrixS=new Matrix() 
GrCMatrixS=new Matrix() 
MFMatrixS=new Matrix() 

objref sav_MFMatrixS,sav_GrCMatrixS,sav_GoCMatrixS
sav_MFMatrixS=new File()
sav_GrCMatrixS=new File()
sav_GoCMatrixS=new File()

//Open recording Files
sav_MFMatrixS.wopen("MFMatrixS.dat")
sav_GrCMatrixS.wopen("GrCMatrixS.dat")
sav_GoCMatrixS.wopen("GoCMatrixS.dat")


//Data for Raster Plots
//GoC
objref timevec1, idvec1, recncs1, tobj1, nil1
timevec1 = new Vector()
idvec1 = new Vector()
recncs1 = new List()
for i=0, NumberCellsd-1 {
  GoC[i].soma tobj1 = new NetCon(&v(0.5), nil1)
  tobj1.record(timevec1, idvec1, i+1) 
  recncs1.append(tobj1)
}
objref tobj1 

//GrC
objref timevec2, idvec2, recncs2, tobj2, nil2
timevec2 = new Vector()
idvec2 = new Vector()
recncs2 = new List()
for i=0,NumberCellsb-1 {
  GrC[i].soma tobj2 = new NetCon(&v(0.5), nil2)
  tobj2.record(timevec2, idvec2, i+1) 
  recncs2.append(tobj2)
}
objref tobj2