//
//  Read spike data from spiking neuron network simulations
//  - draw rasterplot for selected cells
//  - distribute cells randomly on a 4mm x 4mm plane
//  - calculate LFP with the kernel method
//
//  Written by Alain Destexhe, CNRS, 2020
//

load_file("nrngui.hoc")


// 1. read data in vectors

N = 5000	// nb of cells to consider
Ne = 4000	// nb of excitatory cells
Ni = 1000	// nb of inhibitory cells

NSe = 10133	// nb of excitatory spikes
NSi = 9728	// nb of inhibitory spikes
totspikes = NSe+NSi	// total nb of spikes
tmin = 9000	// min time (to skip)
tmax = 1000	// max time

objref Spikes[N]	// vectors to store spikes
for i=0,N-1 {
  Spikes[i] = new Vector(2500)
}

double nspikes[N]	// vector to store the nb of spike of each neuron

ropen("brunel_exc.txt")	// open file for excitatory cells

for i=0,NSe-1 {
  ncell = fscan()               // read next cell spiking
  time = 1000 * fscan()-tmin		// read next spike time
  if(ncell<N) {
    j = nspikes[ncell]
    Spikes[ncell].set(j,time)	// store spike time
    nspikes[ncell] = j+1
  }
}


ropen("brunel_inh.txt")	// open file for inhibitory cells

for i=0,NSi-1 {
  ncell = Ne + fscan()		// read next cell spiking
  time = 1000 * fscan()-tmin	// read next spike time
  if(ncell<N) {
    j = nspikes[ncell]
    Spikes[ncell].set(j,time)	// store spike time
    nspikes[ncell] = j+1
  }
}


// 2. draw raster

Nmin = 0		// first cell to draw in rasterplot
Nmax = N-1		// last cell to draw in rasterplot
Nstp = 10		// step cell to draw


objectvar g1		// create graph
g1 = new Graph()
g1.size(0,tmax,0,N+1)

for (i=Nmin;i<=Nmax;i=i+Nstp) {	// loop on cells
  for j=0,nspikes[i]-1 {  	// loop on spikes
    time = Spikes[i].get(j) 	// get spike time
    if(i<Ne) {
	g1.mark(time,i+1,"O",3,3,1)	// draw exc spike
    } else {
	g1.mark(time,i+1,"O",3,2,1)	// draw inh spike
    }
  }
}

g1.flush()



// 3. distribute cells in a 2D grid

xmax = 0.2	// size of the array (in mm)
ymax = 0.2

double x[N],y[N]  // create vectors for coordinates

objref rnd	  // create random generator
rnd = new Random()

for i=0,N-1 {	  // calculate coordinates to distribute neurons
  x[i] = rnd.uniform(0,xmax)
  y[i] = rnd.uniform(0,ymax)
  // if(i<10) print "cell ",i," , x,y=",x[i],y[i]
}



// 4. calculate LFP
//
// Table of respective amplitudes:
// Layer   amp_i    amp_e
// deep    -2       -1.6
// soma    30       4.8
// sup     -12      2.4
// surf    3        -0.8
//

dt = 0.1	// time resolution
npts = tmax/dt	// nb points in LFP vector
objref LFP
LFP = new Vector(npts)	// create vector for LFP

xe = xmax/2
ye = ymax/2	// coordinates of electrode

va = 200	// axonal velocity (mm/sec)
lambda = 0.2	// space constant (mm)
dur = 100	// total duration of LFP waveform
nlfp = dur/dt	// nb of LFP pts
amp_e = 0.7	// uLFP amplitude for exc cells
amp_i = -3.4	// uLFP amplitude for inh cells
sig_i = 2.1	// std-dev of ihibition (in ms)
sig_e = 1.5 * sig_i  // std-dev for excitation

//amp_e = -0.16	// exc uLFP amplitude (deep layer)
//amp_i = -0.2	// inh uLFP amplitude (deep layer)

amp_e = 0.48	// exc uLFP amplitude (soma layer)
amp_i = 3	// inh uLFP amplitude (soma layer)

//amp_e = 0.24	// exc uLFP amplitude (superficial layer)
//amp_i = -1.2	// inh uLFP amplitude (superficial layer)

//amp_e = -0.08	// exc uLFP amplitude (surface)
//amp_i = 0.3	// inh uLFP amplitude (surface)



// calculate the contribution of excitatory cells

s_e = 2*sig_e*sig_e
s_i = 2*sig_i*sig_i

i=0
for i=0,Ne-1 {	// loop on excitatory cells
  dist = sqrt( (x[i]-xe)^2 + (y[i]-ye)^2 )	// calc distance to electrode (in mm)
  delay = 10.4 + dist/va			// delay to peak (in seconds)
  amp = amp_e * exp(-dist/lambda)		// calc LFP amplitude at electrode position
  print "cell ",i," , amplitude = ",amp
  for j=0,nspikes[i]-1 {	// loop on spikes
    time = Spikes[i].get(j)	// get spike time
    tp = time + delay		// peak time of the uLFP
    for k=0,nlfp-1 {		// scan uLFP time
	t = time + k*dt		// time in ms
	kernel = amp * exp(-(t-tp)*(t-tp)/s_e)  // calc kernel
	index = int(t/dt)
	if(index<npts) LFP.set(index, LFP.get(index) + kernel)		// add to vector
    }
  } 
}



// calculate the contribution of inhibitory cells

for u=0,Ni-1 {	// loop on excitatory cells
  i=Ne+u	// index of the cell
  dist = sqrt( (x[i]-xe)^2 + (y[i]-ye)^2 )	// calc distance to electrode (in mm)
  delay = 10.4 + dist/va			// delay to peak (in seconds)
  amp = amp_i * exp(-dist/lambda)		// calc LFP amplitude at electrode position
  print "cell ",i," , amplitude = ",amp
  for j=0,nspikes[i]-1 {	// loop on spikes
    time = Spikes[i].get(j)	// get spike time
    tp = time + delay		// peak time of the uLFP
    for k=0,nlfp-1 {		// scan uLFP time
	t = time + k*dt		// time in ms
	kernel = amp * exp(-(t-tp)*(t-tp)/s_i)  // calc kernel
	index = int(t/dt)
	if(index<npts) LFP.set(index, LFP.get(index) + kernel)		// add to vector
    }
  } 
}





objectvar g2            // create graph
g2 = new Graph()
g2.size(0,tmax,LFP.min(),LFP.max())

LFP.plot(g2,dt)		// plot the LFP