// This models shows how a sequence of 5 pyramidal cells may be replayed in FORWARD and REVERSE direction during sharp-wave ripples.
// A network of 9 clusters, each of 9 pyramidal cells (PC) and 1 basket cell (IN). The network demonstrates a replay of 5 PCs,
// which are bidirectionally synaptically connected in a sequence: 0->3->30->33->60, 60->33->30->3->0 (cell indices).
// The first "cue" afferent EPSP applied to PC(0) initiates the replay sequence, the second cue afferent EPSP is applied to PC(60)
// All 81 PCs are randomly connected into an axonal plexus by gap junctions in their proximal collaterals. 
// The collaterals are stimulated at random, once per stim_interval_percell=20 ms on average.
// Please be patient, simulation may take a few hours, because of the fine spatial discretization of the axons.
// author: Nikita Vladimirov, PhD <nikita.vladimirov@gmail.com>

load_file("nrngui.hoc")
load_file("template_PC.hoc")
load_file("template_IN.hoc")
load_file("gap_collat.hoc")
// size of clusters:
W_neighborhood = 3
H_neighborhood = 3
// overall size of the network of pyramidal cells
Nx=W_neighborhood*3
Ny=H_neighborhood*3
ncells=Nx*Ny // pyramidal cells total
Nbaskets = 9 // basket cells total
// parameters of axonal plexus:
gj_per_cell = 3
gj_radius_x = gj_radius_y = Nx
gj_pos_min = 0.25
gj_pos_max = 0.75
gj_conductance = 3e-3 //[mcS]
stim_interval_percell = 20 // [ms], mean interstimulus interval, per cell
// random number seeds:
stim_seed = 1 // seed for axonal stimulation times
gj_radius_seed = 2 // seed for connecting gap junctions
gj_pos_seed = 3 // seed for gap junctions on the collateral
// synaptic weights:
wsyn_2PC = 0.005 //CA3->pyr
wsyn_PC2PC = 0.0045 //pyr->pyr
wsyn_PC2IN = 0.0035 //pyr->bask
wsyn_IN2PC = 0.070 //bask->pyr; 0.050 default
objref cells, baskets, gaps // cells are pyramidal cells; baskets are the interneurons;
objref PC2PCsynL, PC2INsynL, P2BaskConL, IN2PCsynL, NetConL, NetStimL
objref Ipulse3L, RandNetStimL, CA3NetStim
NetConL = new List()
PC2PCsynL = new List()
PC2INsynL = new List()
P2BaskConL = new List()
IN2PCsynL = new List()
NetStimL = new List()
Ipulse3L = new List()
RandNetStimL = new List()
// coordinates of the extracellular electrode (um):
Ex = (Nx/2)*40 
Ey = 20
Ez = (Ny/2)*40
//////////////////// procedures and functions //////////////////////	
proc mkcells(){ local i localobj pyr
// $1 is the total number of cells
   cells=new List()
   for i=0,$1-1{
    pyr = new pyramidalCA1()
	pyr.geom_nseg() // this is very important after changing geometry
	forsec pyr.axonal { nseg *= 9 }
	pyr.setcurrentbias_soma(-0.05)
	pyr.setcurrentbias_AIS(-0.1)
	cells.append(pyr)
   }
}

proc mkbaskets(){ local i localobj cell_
// $1 is the total number of basket cells
   baskets = new List()
   for i=0,$1-1{
     cell_ = new basket()
	 baskets.append(cell_)
   }
}

proc setpositions(){ local i, x, y, x_index, y_index
  for i=0,cells.count()-1{
	  x_index =i%Nx+1
	  x = (x_index-1)*40
	  y_index =(i+1-x_index)/Nx+1
	  y = (y_index-1)*40
	  cells.object(i).position(x,0,y)
	  cells.object(i).setindexXY(x_index,y_index)
  }
  for i=0,baskets.count()-1 {
	x = (i%3)*40*3 + 40/2
	y = (i-i%3)*40*3 +40/2
	baskets.o(i).position(x,0,y) 
  }
  
}

proc rotateRandOY(){ local i, choices localobj rand
	choices = 4
	rand = new Random()
	rand.discunif( 0, choices - 1 )
	for i=0,cells.count()-1 {
		cells.o(i).rotateOY(2*PI*rand.repick()/choices)
	}
}

obfunc shuffle(){ local n, rand_ind, a localobj rand, shuff
  rand=new Random()
  shuff=new Vector($o1.size(),0)
  for (n=0; n<=$o1.size()-1;n=n+1){
   shuff.x[n]=$o1.x[n]
  }
  for (n=$o1.size(); n>1; n=n-1){
   rand_ind=rand.discunif(0,n-1)
   a=shuff.x[n-1]
   shuff.x[n-1]=shuff.x[rand_ind]
   shuff.x[rand_ind]=a
   }
  return shuff
}

func round() { return int($1+0.5) }

/////////////// Connect the cells into a network by gap junctions ///////////////
proc connect_gj(){ local i, i_gj, icell, xi, yi, xj, yj, n1, n2, attempt, check1, check2, pos1, pos2 localobj cell_gids, cell_gids_shuffled, gap_, stubs, rand2, rand3, gj_table
ngaps=int(ncells*gj_per_cell/2)
stubs=new Vector(ncells+1)// Fortan-style indexing, 1..ncells
gj_table=new Matrix(ngaps+1,2)
rand2=new Random(gj_radius_seed)
rand3=new Random(gj_pos_seed)
rand2.discunif(-gj_radius_x,gj_radius_x)
rand3.uniform(0,1)
stubs.fill(gj_per_cell) //regular
stubs.x[0]=0 // fortran style

gaps=new List()
cell_gids=new Vector(cells.count())
cell_gids.indgen(1,cells.count(),1)
cell_gids_shuffled=shuffle(cell_gids)

i_gj=0
for i=0,cells.count()-1{
  icell=cell_gids_shuffled.x[i]// icell=1..ncells, Fortran style
  xi=cells.object(icell-1).x_index
  yi=cells.object(icell-1).y_index
  while(stubs.x[icell]>0) {
       xj=xi+rand2.repick()
       yj=yi+rand2.repick()
       n1=icell
       n2=Nx*(yj-1)+xj
       attempt=0
	   check1=((xj<1)||(yj<1)||(xj>Nx)||(yj>Ny)||((xj==xi)&&(yj==yi)))
	   if(check1==0){ //if cell is within range, try its free stubs
			check2=(stubs.x[n2]==0)
		}
       while ((check1 || check2) && (attempt<1000)){
	       attempt=attempt+1
           xj=xi+rand2.repick()
           yj=yi+rand2.repick() 
		   n2=Nx*(yj-1)+xj	
	       check1=((xj<1)||(yj<1)||(xj>Nx)||(yj>Ny)||((xj==xi)&&(yj==yi)))
	       if(check1==0){ //if cell is within range, try its free stubs
	       check2=(stubs.x[n2]==0)
	     }
	   }
	   if(attempt<1000 && i_gj<ngaps) {
 	     i_gj=i_gj+1
         gj_table.x[i_gj][0]=n1
         gj_table.x[i_gj][1]=n2
		 pos1 = rand3.uniform(gj_pos_min,gj_pos_max)
		 pos2 = rand3.uniform(gj_pos_min,gj_pos_max)
		 gap_=new gapjunction(cells.object(n1-1),0,cells.object(n2-1),0,gj_conductance,pos1,pos2)
		 gaps.append(gap_)
		 stubs.x[n1]=stubs.x[n1]-1
         stubs.x[n2]=stubs.x[n2]-1
	   } else {
	     stubs.x[icell]=0
	   }
    }
  }
}

func quotent() {
//finds quotent q of a divided by b:
// a = b*q + rem
// $1 is a, $2 is b
return ($1-$1%$2)/$2
}

func ipyramid() { local icol, jrow
//finds the first i of pyramidal cell in cluster of 9 cells
// $1 is the index of basket cell
icol = ($1%W_neighborhood)*W_neighborhood
jrow = quotent($1,H_neighborhood)*H_neighborhood
return jrow*Nx + icol
}

func ipyramids() {
// returns the absolute indices of pyr. neighbors of $1 = 0..8 basket cell, 
// by using relative index of pyr. cell in neighborhood $2 = 0..8
return ipyramid($1) + $2%3 + quotent($2,W_neighborhood)*Nx
}

proc connectPCs_to_IN() { local i,j,j0 localobj syn_, netcon_
	for i = 0, baskets.count()-1 {
			baskets.o(i).soma syn_ = new AlphaSyn(0.5)
			syn_.tau = 0.8 // 0.8, Traub et al 2005, AMPA pyr->FS
			for j0 = 0, 8 { //run thru all pyrs in the neighborhood
				j = ipyramids(i,j0) 
				cells.o(j).collat[1] netcon_ = new NetCon(&v(1),syn_)
				netcon_.weight = wsyn_PC2IN
				netcon_.delay = 1
				P2BaskConL.append(netcon_)
				PC2INsynL.append(syn_)
			}
			}
}

proc connectIN_to_PCs() { local j localobj syn_, netcon_
		for i = 0, baskets.count()-1 {
			for j0 = 0, 8 { //run thru all pyrs in the neighborhood
			j = ipyramids(i,j0) //retrieve pyr's gloal index
			cells.o(j).soma syn_ = new Exp2Syn(0.5)
			syn_.tau1 = 2
			syn_.tau2 = 5
			syn_.e = -75
			baskets.o(i).soma netcon_ = new NetCon(&v(1),syn_)
			netcon_.weight = wsyn_IN2PC
			netcon_.delay = 1
			NetConL.append(netcon_)
			IN2PCsynL.append(syn_)
			}
		}
}

proc connectPC_to_PC() { local i localobj syn_, netcon_
// $1 is pre-synaptic, $2 is post-synaptic cell
		 cells.o($2).basal[4] syn_=new AlphaSyn(0.5)
		 syn_.tau = 2 //ms
		 cells.o($1).collat[1] netcon_ = new NetCon(&v(1),syn_)
		 netcon_.weight = $3
		 netcon_.delay = 1
		 PC2PCsynL.append(syn_)
		 NetConL.append(netcon_)
}

proc cueEPSP() { local i localobj syn, netcon, netstim
		netstim = new NetStim()
		netstim.number = 1 // 
		netstim.start = $2
		cells.o($1).shaft syn = new AlphaSyn(0.5)
		syn.tau = 2
		netcon = new NetCon(netstim,syn)
		netcon.weight = wsyn_2PC
		netcon.delay = 1
		NetStimL.append(netstim)
		PC2PCsynL.append(syn)
		NetConL.append(netcon)
}

proc stim_all_collat_random() { local i localobj pulse_, netstim_, netcon_
 for i = 0, cells.count()-1 {
	cells.o(i).collat {
		pulse_ = new Ipulse3(0.9)
		pulse_.dur = 1
		pulse_.amp = 0.05
		netstim_ = new NetStim()
		netstim_.interval = stim_interval_percell
		netstim_.number = 1e9 // essentially unlimited
		netstim_.start = 5
		netstim_.noise = 1 // 0 for deterministic, 1 for negexp
		netstim_.seed(i+stim_seed)
		netcon_ = new NetCon(netstim_,pulse_)
		netcon_.threshold = 0
		netcon_.weight = 1
		netcon_.delay = 0
		Ipulse3L.append(pulse_)
		RandNetStimL.append(netstim_)
		NetConL.append(netcon_)
		}
 }
}


/////////////////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////   execution   /////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////////////////

mkcells(Nx*Ny)
mkbaskets(Nbaskets)
setpositions()
rotateRandOY()
connect_gj()
connectPCs_to_IN()
connectIN_to_PCs()
cueEPSP(0,20) //$1 is cell index, $2 is time of onset
cueEPSP(60, 150) 
// connect for forward replay:
connectPC_to_PC(0,3,wsyn_PC2PC) 
connectPC_to_PC(3,30,wsyn_PC2PC) 
connectPC_to_PC(30,33,wsyn_PC2PC) 
connectPC_to_PC(33,60,wsyn_PC2PC) 
// connect for reverse replay:
connectPC_to_PC(3,0,wsyn_PC2PC) 
connectPC_to_PC(30,3,wsyn_PC2PC) 
connectPC_to_PC(33,30,wsyn_PC2PC) 
connectPC_to_PC(60,33,wsyn_PC2PC) 
//
stim_all_collat_random()

load_file("LFP.hoc")
load_file("81PC_9INforw_rev.ses")