Synch_DELAY_MdB.hoc

load_file("nrngui.hoc")
load_file("cellM1.hoc")
load_file("correl.hoc")
cvode_active(1)

//definiamo la temperatura di funzionamento, 
//altrimenti neuron userebbe sempre quella del calamaro giagante pari a 6 gradi
celsius=34

// Migliore-Hines synchronization and gap junctions. 
// olfactory bulb _correlation 
// Tom & mc Tavish: sync cambiando la posizione delle sinapsi
// i_membrane_ in CVode.use... insieme delle correnti in ingresso nel soma

// migliore ascoli 2005 fixnseg (per le divisioni )

// CIRCUITO TRISINAPTICO
// 20 40 60 80 100Hz Poisson basali 5 morfologie ottimizzate

//definizione delle variabili

// relazione tra neuroni inibitori ed eccitatori. sul soma
// Megias 2002 conta i numeri delle sinapsi
// quanto interneuroni stanno contribuendo sul soma
// grossezza dei dentriti 1.2 um

objref stim[2], syn[2], apc[2]
objref synD[2][400], nt[2][2], nsD[2][400], ncD[2][400], ns, ncS[2], randDen
objref fseq, fv, finp,fout,fphi,fspk,fcorr,fc

objref b, c, d, e
objref b2,c2,d2,e2
// stim, distrx, distry, cdistry, p, rand,syn[50],ns[50],nci[50]
// Esti = new Vector(2,0) // extra neuron spike times interval vettori di 2 elementi inizializzati in 0

double  vecAD[2][400], DistSD[2][400], DistSDm[2], dSDm[2], ds[2]
double  tspk[2][1000], nspk[2], nspkm[2], f[2], Esti[1000], Esi[1000], corr[10000000], corrvi[10000000], tt[10000000]
double  tspkprec[2],tspknext[2],i_spk[2], phi[2]
strdef llab
strdef fphin,foutn,fvn,fcn,ftn,fcrn

objref C1, C2, tC1, tC2

  Nexp    = 1
  wi      = 0.05
  ws      = 5e-4
  Nsyn    = 20
  ds[0]   = 100
  ds[1]   = 150
  Dsd     = 25
  chnoise = 1
  FreqInp = 60
  del     = 2

print "Nsyn=", Nsyn," wi=",wi," ws=",ws," ds0=",ds0," ds1=",ds1," Dsd=",Dsd," Noise=", chnoise," f=", FreqInp," Del=", del

SynLoc = 0.5
tstop=1000
cvode.event(tstop)

sprint(fcn,"Synchro_C_wi%g_dl%g.out",wi,del)
fc = new File(fcn)
fc.wopen()
fc.printf("# Sm     Cm      C0      Dphim\n")

sprint(ftn,"Synchro_Tspk_wi%g_dl%g.out",wi,del)
fspk = new File(ftn)
fspk.wopen()
fspk.printf("# Ne  Nsyn  ws  wi    nspk[0] tspk[0]   nspk[1] tspk[1]\n")

xseed=16327
print xseed

N   = 2    //Num neurons
objref NC[2]
for j=0, N-1{
	NC[j] = new CA1_PC_cAC_sig() 
	NC[j].init()
}

Napic = 189 

for j=1, N-1 {
	NC[j].soma {
	    print secname(), " Nsecs: ", n3d()
    	for i=0,n3d()-1 { 
	    	pt3dchange(i, x3d(i)+400, y3d(i), z3d(i), diam3d(i))
		}
	}
}

ns = new NetStim(0.5)      //Unique current pulse for all the synapses
ns.number = 1000000
//ns.interval = 20          //20ms -> 1/(20e-3) = 50 Hz
ns.interval = 1000/FreqInp  //interval in ms
ns.start = 10
ns.noise = chnoise

for m = 0, (N-1) NC[m].soma {
    syn[m] = new ExpSyn(0.5)
	syn[m].tau=30
	syn[m].e  =-80
	
    nspk[m] =0
    nspkm[m]=0

    apc[m] = new APCount(.5)
    apc[m].thresh=0
}

for m = 0, (N-1) {
	for k = 0,(N-1) {
		NC[m].soma nt[m][k] = new NetCon(&v(.5), syn[k], 0, del, wi)
	}
}

mcell_ran4_init(xseed)
rrr = mcell_ran4(&xseed)
print rrr

for l = 0, N-1 {
   	tspkprec[l] = 0
	tspknext[l] = 0
    i_spk[l]= 0
	l2=1
   	while ( tspk[l][l2] > 0) {
	    tspk[l][l2] = 0
	    l2 = l2 + 1
	} 
}

for i = 0,999 {
    Esti[i]=0.
    for j=0,N-1 {
	   tspk[j][i]=0.
    }
}

for i = 0,399 {
    for j=0,N-1 {
	   vecAD[j][i] = 0
	   DistSD[j][i] = 0
    }
}


kcommon=0
mincommon = 0
NspikesTot =0
deltats = 0.0

b = new VBox()
b.intercept(1)

	c = new Graph()
	c.size(0,2000,-1,1)
	c.xaxis(1)
	c.exec_menu("10% Zoom out")
	c.color(2)
	c.label(0.43,0.93,"Corr V1-V2")
	c.label(0.49,0.0,"t")

	d = new Graph()
	d.size(0,2000,0,1)
	d.xaxis(1)
	d.exec_menu("10% Zoom out")
	d.color(3)
	d.label(0.43,0.93,"S(t)")
	d.label(0.49,0.0,"t")
		
	e = new Graph()
	e.size(0,2000,0,50)
	e.xaxis(1)
	e.exec_menu("10% Zoom out")
	e.color(5)
	e.label(0.45,0.93,"N_spikes")
	e.label(0.49,0.0,"t")
	
	xpanel("")
	xbutton("Shuffle syn. & reRun ", "runm()")
	xpanel()
	xpanel("variables")
	xvalue("inter-soma Synaptic Inhibition - wi","wi",0,"wca()")
	xvalue("Dendr Net Synaptic Conduptance - ws","ws",0,"wca()")
	xpanel()
	
b.intercept(0)
b.map()

b2 = new VBox()
b2.intercept(1)

	c2 = new Graph()
	c2.size(0,1,-1,1)
	c2.xaxis(1)
	c2.exec_menu("10% Zoom out")
	c2.color(2)
	c2.label(0.4,0.93,"<Corr V1-V2>")
	c2.label(0.49,0.0,llab)
		
	d2 = new Graph()
	d2.size(0,1,0,1)
	d2.xaxis(1)
	d2.exec_menu("10% Zoom out")
	d2.color(3)
	d2.label(0.4,0.93,"<S(t)>")
	d2.label(0.49,0.0,llab)
		
	e2 = new Graph() 
	e2.size(0,1,0,55)
	e2.xaxis(1)
	e2.exec_menu("5% Zoom out")
	e2.color(5)
	e2.label(0.42,0.93,"<Freq.>")
	e2.label(0.49,0.0,llab)
	
b2.intercept(0)
b2.map()



proc advance() {    //Neuron time process
	fadvance()      //each dt increment 

    for l = 0, N-1 {
		if ( apc[l].n > nspk[l] ) {  
            
    		 nspk[l] = apc[l].n 
    		 tspk[l][nspk[l]] = t
    		 NspikesTot = NspikesTot + 1

             if (l == 0) fspk.printf("%d %d %g %g   %d %g   %d %g\n",Ne,Nsyn,ws,wi,nspk[0],t,nspk[1],0)
             if (l == 1) fspk.printf("%d %d %g %g   %d %g   %d %g\n",Ne,Nsyn,ws,wi,nspk[0],0,nspk[1],t)

    		 if (l == 0) e.mark(t,nspk[l],"o",8,1,1)
    		 if (l == 1) e.mark(t,nspk[l],"t",8,1,1)
    		 e.flush
    	}
	}

    istep = istep+1
    sum12 = sum12 + NC[0].soma.v*NC[1].soma.v*dt
    sum1  = sum1  + NC[0].soma.v^2.*dt
	sum2  = sum2  + NC[1].soma.v^2.*dt
    
    if ( istep == NstepsCorr ){
       cc   = sum12/sqrt(abs(sum1*sum2))
//       print Nwindowscorr, "cc ", cc, "ccvi ",ccvi, "sumI ",sumI, "sum1 ",sum1, "sumVI ",sumVI, "s_m ",s_med //, "c_m ",corr_med
       corr[ic] = cc
       tt[ic] = t
       
       c.color(1)
       c.line(ic,cc)
       c.flush
              
       ic = ic + 1
       istep = 0
       sum1  = 0
       sum2  = 0
       sum12 = 0
       Nwindowscorr = Nwindowscorr + 1
       corr_med   = corr_med   + cc
    } 
    if ( Ne == 1 ) {
        fv.printf("%g %g %g %g \n", t, NC[0].soma.v, NC[1].soma.v, NC[1].apic[indAD].v)
    }

    tp=t
}


proc wca() {
	for k = 0, N-1 {
		for m = 0, N-1 {
      		nt[k][m].weight=wi
		}
	}
}

proc CreateShuffleSyn() {

for m=0,N-1 {
	access NC[m].soma
	distance()
	dm=0
	for j=0,(Nsyn-1)  {

		k=0
		while(k==0){
    		indAD = int(Napic*mcell_ran4(&xseed)) 
 			rrr = mcell_ran4(&xseed)

 			NC[m].apic[indAD] {
 				if ( distance(rrr) > ds[m] && distance(rrr) < (ds[m]+Dsd) && diam(rrr) < 1.2 ) {
					k=1
//  				print indAD, distance(rrr) 
    	    	}
			}
		}
	    access NC[m].apic[indAD] // random dendrites

	   	synD[m][j] = new Exp2Syn(rrr)     // nuova sinapsi a due tempi
   		synD[m][j].tau1 = 0.5             // costante di attivazione   (ms)
	   	synD[m][j].tau2 = 3               // costante di inattivazione (ms)
   		synD[m][j].e    = 0               // sinapsi di tipo eccitatorio
   		ncD[m][j] = new NetCon(ns,synD[m][j],0,0,ws) //Unico impulso di corrente per tutte le sinapsi

//		print m, j, indAD, synD[m][j].get_loc(),rrr

	    DistSD[m][j] = distance(rrr)
		dm = dm + DistSD[m][j]
	}
	DistSDm[m] = dm/Nsyn
}
print "      Synapsis shuffling - dSm[0] = ",DistSDm[0], " - dSm[1] = ",DistSDm[1] 
print "                   Frequency (Hz) = ",1000/ns.interval, " - ISI (ms) = ",ns.interval
print "                               ws = ",ws, " - wi = ",wi, " - delay = ",del
}


proc runMain() {

    CreateShuffleSyn()
	
    Corr_me   = 0.0
    CS_me      = 0.0
    CS2_me     = 0.0

	ic  = 0
	cc  = 0

	for k = 0, N-1 {

		for j=0,Nsyn-1 {
	    	ncD[k][j].weight = ws
		}
		for m = 0, N-1 {
      		nt[k][m].weight=wi
		}
		apc[k].n      = 0
		f[k]          = 0.0
		nspkm[k]      = 0
 	}

    sprint(fphin,"Synchro_Phi_wi%g_dl%g.out",wi,del)
	fphi = new File(fphin)
	fphi.wopen()
	fphi.printf("# Ne,Nsyn,ws,wi,nspk[0],i_spk[0],tspkprec[0],tspknext[0],phi[0],    nspk[1],i_spk[1],tspkprec[1],tspknext[1],phi[1],    t, s_med, cs_med \n")

    sprint(fvn,"Synchro_V_wi%g_dl%g.out",wi,del)
	fv = new File(fvn)
	fv.wopen()
	fv.printf("# t    NC[0].soma.v, NC[1].soma.v  NC[1].apic[%d].v \n",indAD)

	if ( Ne == 1 ) {
		fv.printf("# Nsyn %d  - d0 %g  - d1 %g  - ws %g  - wi %g  - dli %g \n", Nsyn,DistSDm[0],DistSDm[1],ws,wi,del)
		fc.printf("# Nsyn %d  - d0 %g  - d1 %g  - ws %g  - wi %g  - dli %g \n", Nsyn,DistSDm[0],DistSDm[1],ws,wi,del)
    }
    print "kakbNeNsy N0 N1 - ws  wi  f1  f2 Crm Sm dPhi0 dPhi1 CSm"

   	for k = 0, N-1 {
		dSDm[k] = 0
	}
    DPhim  = 0.0
    DPhim2 = 0.0

	for ke = 0, Nexp -1 {

        Ne = ke + 1 
    	Nwindowscorr=0
	    istep = 0
    	sum1  = 0.0
	    sum2  = 0.0
    	sum12 = 0.0
		
		if( ke > 0 ) {
	    	for k = 0, N-1 {
				DistSDm[k] = 0
			}
			CreateShuffleSyn()
		}
    	
		run()
		
		for k = 0, N-1 {
            dSDm[k] = dSDm[k] + DistSDm[k]
			if ( apc[k].time > 0 ) {f[k] = f[k] + nspk[k]/apc[k].time*1000}
			nspkm[k] = nspkm[k] + nspk[k]
			NspikesTotm = NspikesTotm + nspk[k]
		}
		
		d.beginline()
        tmax = tspk[0][nspk[0]]
        if ( tmax > tspk[1][nspk[1]] ) { tmax = tspk[1][nspk[1]] }
		ttt = 0
		ddt = 0.1
		Ntsteps = 0
		DeltaPhim = 0
		
    	for l = 0, N-1 {
			i_spk[l] = 0
			tspknext[l]  = tspk[l][i_spk[l]+1]
		}
		
   		phi[0]=mcell_ran4(&xseed)*6.283184   
   		phi[1]=mcell_ran4(&xseed)*6.283184   
		cs_med = 0.0

		while ( ttt < tmax ) {

	   		phi[0]=mcell_ran4(&xseed)*6.283184
   			phi[1]=mcell_ran4(&xseed)*6.283184

			if (nspk[0] >1 && nspk[1] >1) {
	    	for l = 0, N-1 {

		        if ( ttt >= tspk[l][i_spk[l]+1] ) {
		            i_spk[l] = i_spk[l] + 1
		            tspkprec[l] = tspk[l][i_spk[l]] 
		            tspknext[l] = tspk[l][i_spk[l]+1] 
	    		}

		  	    if ( (tspkprec[l] != tspknext[l]) ) {
		    		phi[l]=0    // rnd:22/05/2023
		  	        if  (ttt < tspknext[l]) {
			  	        phi[l] = 2.*3.141592 * (ttt - tspkprec[l]) / (tspknext[l] - tspkprec[l])
					}
    			}
    		}
			}
		
    		// in each step, updates the phase difference
    		DeltaPhim = DeltaPhim + (phi[1]-phi[0])*ddt
			cs_med =  cs_med + cos( (phi[1]-phi[0]) ) * ddt
		    ttt = ttt + ddt
		    Ntsteps = Ntsteps + 1
//     	    print "tsp[0]=", tspkprec[0], " tsn[0]=", tspknext[0]," phi[0]=", phi[0],"tsp[1]=", tspkprec[1], " tsn[1]=", tspknext[1]," phi[1]=", phi[1], "s(t)=", s_med/ttt
            
            if ( Ne == 1 ) {  
	            if ( Ntsteps % 10 == 0 ){
    	    		fphi.printf("%d %d %g %g %d %d %g %g %g    %d %d %g %g %g   %g %g   %g %g   %g %g %g\n",Ne,Nsyn,ws,wi,nspk[0],i_spk[0],tspkprec[0],tspknext[0],phi[0],nspk[1],i_spk[1],tspkprec[1],tspknext[1],phi[1], ttt,s_med/ttt,dSDm[0],dSDm[1],DeltaPhim/ttt,cs_med/ttt,Corr_me)
				}
				d.line(ttt,s_med/ttt)
    		    d.flush
            }

		}
		
//	   		print tmax, nspk[0], nspk[1]
        if (tmax == 0) {
            ttt=1
	        cs_med = cos( (phi[1]-phi[0]) )
	        DeltaPhim = (phi[1]-phi[0])
        }
		CS_me   = CS_me + cs_med/ttt
		DPhim   = DPhim  + DeltaPhim/ttt
		CS2_me    = CS2_me + cs_med/ttt * cs_med/ttt
		DPhim2    = DPhim2  + DeltaPhim/ttt * DeltaPhim/ttt

        varS  = abs(S2_me/Ne - S_me/Ne * S_me/Ne)
        varCS = abs(CS2_me/Ne - CS_me/Ne * CS_me/Ne)
        varPH = abs(DPhim2/Ne - DPhim/Ne * DPhim/Ne)
		
		S_sig     = sqrt(varS) 
		CS_sig    = sqrt(varCS)
		DPhim_sig = sqrt(varPH)

        fphi.close
        
    	print Ne,Nsyn,nspkm[0]/Ne,nspkm[1]/Ne,ws,wi,f[0]/Ne,f[1]/Ne,Corr_me/Ne,S_me/Ne,dSDm[0]/Ne,dSDm[1]/Ne,DPhim/Ne,CS_me/Ne,S_sig,CS_sig,DPhim_sig,del
    	print "           Sm = ", S_me/Ne,"  Cm = ", CS_me/Ne,"  C(0) = ", Corr_me/Ne

         
	    NspikesTot= 0
    	for l = 0, N-1 {
 
    		nspk[l]=0
	    	tspkprec[l] =0
    		tspknext[l] =0
    		l2=1
	    	while ( tspk[l][l2] > 0) {
			    tspk[l][l2] = 0
			    l2 = l2 + 1
			}
			apc[l].n=0 
    	}

    	ic = 0

	    sprint(foutn,"Synchro_wi%g_dl%g.out",wi,del)
		fout = new File(foutn)
		fout.wopen()
		fout.printf("#  ds1  NspkT  nspk[0] nspk[1]  Nsyn  ws   wi   f1m   f2m   corr_vv   S_me   dS0  dS1  Dphinm   CS_me    S_sig  CS_sig   DPhim_sig  Nexp del\n")
		fout.printf("%g %g %g %g %d %g %g %g %g %g %g %g %g %g %g %g %g %g %d %g \n", ds1,NspikesTotm/Ne,nspkm[0]/Ne,nspkm[1]/Ne,Nsyn,ws,wi,f[0]/Ne,f[1]/Ne,Corr_me/Ne,S_me/Ne,dSDm[0]/Ne,dSDm[1]/Ne,DPhim/Ne,CS_me/Ne,S_sig,CS_sig,DPhim_sig,Ne,del)
		fout.close()

        ccc=0
	    if ( Nwindowscorr > 0 ) ccc = corr_med  / Nwindowscorr
	    fc.printf("%g %g %g %g \n",s_med/ttt,cs_med/ttt,ccc,DeltaPhim/ttt)
	    
		if ( Ne == 1 ) {
			fv.printf("# %g %g %g %g %d %g %g %g %g %g %g %g %g %g %g %g %g %g %d %g \n", ds1,NspikesTotm/Ne,nspkm[0]/Ne,nspkm[1]/Ne,Nsyn,ws,wi,f[0]/Ne,f[1]/Ne,Corr_me/Ne,S_me/Ne,dSDm[0]/Ne,dSDm[1]/Ne,DPhim/Ne,CS_me/Ne,S_sig,CS_sig,DPhim_sig,Ne,del)
	    	fv.close()
	   	}
	   	
	}

   	fc.close()

	xval = ds1
   	c2.mark(xval,Corr_me/Nexp,"o",8,1,1.5)
	c2.label(0.49,0.0,llab)
   	c2.flush

   	d2.mark(xval,S_me/Nexp,"o",8,1,1.5)
	d2.label(0.49,0.0,llab)
   	d2.flush


   	e2.mark(xval,f[0]/Nexp,"o",8,1,1.5)
   	e2.mark(xval,f[1]/Nexp,"t",8,1,1.5)
	e2.label(0.49,0.0,llab)
   	e2.flush

	fout.close()
	fspk.close()

}   // fine runMain

load_file("Ses_t1000.ses")

runMain()