Two Models for synaptic input statistics for the MSO neuron model (Jercog et al. 2010)

 Download zip file 
Help downloading and running models
Accession:141991
The model is a point neuron model with ionic currents from Rothman & Mannis (2003) and with an update of the low threshold potassium current (IKLT) measured in-vitro by Mathews & Jercog et al (2010). This model in conjunction with the synaptic input models presented here has been used to gain insight into mechanisms that account for experimentally observed asymmetries in ITD tuning (Brand et al, 2002). Asymmetry and displacement of the ITD response function is achieved in the model by the interplay between asymmetry of the excitatory inputs arriving from the two sides and the precise voltage dependent activation of IKLT. In Jercog et al (2010) we propose two different mathematical ways, physiologically plausible scenarios, of generating the asymmetry in the bilateral synaptic input events. Here, we present two models for simulating the stochastic synaptic input trains.
Reference:
1 . Jercog PE, Svirskis G, Kotak VC, Sanes DH, Rinzel J (2010) Asymmetric excitatory synaptic dynamics underlie interaural time difference processing in the auditory system. PLoS Biol 8:e1000406 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Axon;
Brain Region(s)/Organism:
Cell Type(s): Medial Superior Olive (MSO) cell;
Channel(s):
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: C or C++ program;
Model Concept(s): Synaptic-input statistic; Audition;
Implementer(s): Jercog, Pablo [pjercog at gmail.com];
#include <stdio.h>
#include <math.h>
#include <time.h>
#include <stdlib.h>

//#####################################################################
/* Runge Kutta integrator from numerical recipies plus improvements */
/* void *deri(int n,double h[],double D[],double t);  */
/* function argument not tested yet */

void rk4(void deri(int , double [], double [], double ), \
		 double h[], int n, double t, double dt)
{
#define naux 60 
	
	int i;
	double k1[naux],k2[naux],k3[naux],k4[naux],h0[naux];
	double dt2, dt6;
	
	dt2=dt/2.;
	dt6=dt/6.;
	
	for (i = 0 ; i<n; i++)
		h0[i] = h[i];
	
	deri(n,h0,k1,t);
	for (i =0 ; i<n ; i++)
		h0[i]=h[i]+dt2*k1[i];
	
	deri(n,h0,k2,t+dt2);
	for (i =0 ; i<n ; i++)
		h0[i]=h[i]+dt2*k2[i];
	
	deri(n,h0,k3,t+dt2);
	for (i =0 ; i<n ; i++)
		h0[i]=h[i]+dt*k3[i];
	
	deri(n,h0,k4,t+dt);
	for (i = 0 ; i<n ; i++)
	{h0[i]=h[i]+dt*k4[i];};
	
	for (i =0; i<n ; i++)
		h[i]=h[i]+dt6*(2.*(k2[i]+k3[i])+k1[i]+k4[i]);
	
	return;
}
/*End of Runge Kutta Code*/
//#####################################################################



//#####################################################################
/*Global variables definitions*/


#define Nneurons_exc 128 //256
#define Nneurons_inh 128
#define PI 3.1415926535897932384626433832795028841971 
#define Ntrials 20 //20
#define Npoints 60 //40
#define Nspikes 40 //20 //10
#define N_steps_sta 10000
#define Ndiscret 1000

double V,a,b,c,w,z,nn,p,m,h,r,x,y,xx,yy,xxx,yyy,xxxx,yyyy;
double cur_term,IA,ILT,IHT,INA,IH,ILK,syn_term,syn1ex,syn2ex,syn1in,syn2in,CM;
double a_inf,b_inf,c_inf,w_inf,z_inf,n_inf,p_inf,m_inf,h_inf,r_inf,tau_a,tau_b,tau_c,tau_w,tau_z,tau_n,tau_p,tau_m,tau_h,tau_r;
double tauex1,tauex2,tauinh1,tauinh2,GEXC1,GEXC2,GINH1,GINH2;

void takens(int n, double v[], double dv[], double t)
/* Takens equations */
{
 V=v[0];
 a=v[1];
 b=v[2];
 c=v[3];
 w=v[4];
 z=v[5];
 nn=v[6];
 p=v[7];
 m=v[8];
 h=v[9];
 r=v[10];
 x=v[11];
 y=v[12];
 xx=v[13];
 yy=v[14];
 xxx=v[15];
 yyy=v[16];
 xxxx=v[17];
 yyyy=v[18];
                                                                                            
dv[0]=0.001*(cur_term-IA-ILT-IHT-INA-IH-ILK-syn_term)/CM;
dv[1]=(a_inf-v[1])/tau_a;
dv[2]=(b_inf-v[2])/tau_b;
dv[3]=(c_inf-v[3])/tau_c;
dv[4]=(w_inf-v[4])/tau_w;
dv[5]=(z_inf-v[5])/tau_z;
dv[6]=(n_inf-v[6])/tau_n;
dv[7]=(p_inf-v[7])/tau_p;
dv[8]=(m_inf-v[8])/tau_m;
dv[9]=(h_inf-v[9])/tau_h;
dv[10]=(r_inf-v[10])/tau_r;
dv[11]=v[12];
dv[12]= - 2.*v[12]/tauex1 - v[11]/(tauex1*tauex1);
dv[13]=v[14];
dv[14]= - 2.*v[14]/tauex2 - v[13]/(tauex2*tauex2);
dv[15]=v[16];
dv[16]= - 2.*v[16]/tauinh1 - v[15]/(tauinh1*tauinh1);
dv[17]=v[18];
dv[18]= - 2.*v[18]/tauinh2 - v[17]/(tauinh2*tauinh2);

return;
}

int main(){
FILE *prnt1, *prnt2, *prnt3, *prnt4,*prnt5;
double v[19], e1[Nneurons_exc], e2[Nneurons_exc],e3[Nneurons_inh],e4[Nneurons_inh];
double sta[N_steps_sta],sta_stdev[N_steps_sta],memory[N_steps_sta],iklt_memory[N_steps_sta],gklt_memory[N_steps_sta],epsp_gklt_ave[N_steps_sta],uper[N_steps_sta],lower[N_steps_sta],epsp_ave[N_steps_sta],epsp_ave_lower[N_steps_sta],epsp_ave_upper[N_steps_sta],epsp_stdev[N_steps_sta],iklt_ave[N_steps_sta],gklt_ave[N_steps_sta],ipsp_ave[N_steps_sta];
double syn1[N_steps_sta],syn2[N_steps_sta],syn_ex1[N_steps_sta], syn_ex2[N_steps_sta],syn_spike_ex1[N_steps_sta],syn_spike_ex2[N_steps_sta];
int old1[Nneurons_exc],old2[Nneurons_exc],old3[Nneurons_inh],old4[Nneurons_inh];
double T_exc_ipsi[Nneurons_exc][Nspikes],T_exc_contra[Nneurons_exc][Nspikes],T_inh_ipsi[Nneurons_inh][Nspikes],T_inh_contra[Nneurons_inh][Nspikes];
double A_exc_ipsi[Nneurons_exc][Nspikes],A_exc_contra[Nneurons_exc][Nspikes],A_inh_ipsi[Nneurons_inh][Nspikes],A_inh_contra[Nneurons_inh][Nspikes];
double I0, IEXT,q0p17,q3p03;
double startex,EEXC,startinh,EINH,ZETA,PHI,tauex,GEXC;
double EK,ENA,EH,ELK,GNA,GHT,GLT,GA,GH,GLK,VREST,RREST,TAUM,VTH,SM,GE22,GE38;
double t,dt,tTOT,T,freq,v0,ITDmax;
double u1,u2,R,Phi;
double sigmaT1,sigmaT2,sigmaT3,sigmaT4,sigmaA1,sigmaA2,sigmaA3,sigmaA4,v1,v2,v3,integral;
double release_now,inh_ipsi_delay,exc_contra_delay,inh_contra_delay,exc_ipsi_delay;
double prob_max_ipsi_exc,prob_max_contra_exc,prob_max_ipsi_inh,prob_max_contra_inh,threshold1,threshold2,threshold3,threshold4,lambda1,lambda2,lambda3,lambda4,nu1,nu2,nu3,nu4;
int i,j,k,l,o,h,p,n_cycles;
int N,decision,Npeak,N_over_trials,P;
int Nsmallepsp1,Nsmallepsp2,Nsmallepsp3,Nsmallepsp4,k1,k2,k3,k4;
double R1x,R1y,R2x,R2y,R3x,R3y,R4x,R4y,R1,R2,R3,R4,y2,R_ss;
double Gtotal, gA,gHT,gLT,gNA,gH,gLK;
double V0,w0,z0,n0,p0,m0,h0,r0;


prnt1=fopen("Model1_STA_epsg_epsp_iklt_w-inh.dat","w");//ITD_w-IKLT_09_-099_09_500.dat
fclose(prnt1);
prnt2=fopen("Model1_train_EPSGs-EPSPs_w-inh.dat","wt");
prnt3=fopen("Model1_ITD-curve_w-inh.dat","w");//ITD_w-IKLT_w-INH_09_-099_09_500.dat
fclose(prnt3);

srand48(time(0));

//*********************************//
//Definition os the local constants//
//*********************************//																						
freq=.5;//in KHz

I0=0.0;
IEXT=0.0;
q0p17=0.17;//1.0;
q3p03=3.03;//1.0;
                                                                                            
//parameters for syn_term:
syn1ex=0.;
syn2ex=0.;
syn1in=0.;
syn2in=0.;
tauex1=0.1;
tauex2=0.1;
tauinh1=0.4;
tauinh2=0.4;
GEXC1=2.25;//2.25//0.7;//15.;//0.72;//.55;//.50;//1.35;
GEXC2=2.25;//2.25//(GEXC1*tauex1)/tauex2;
GINH1=0.0;//0.05;//0.05;//.70;
GINH2=1.4;//1.2;//0.40;//0.17;//0.05;//.7;

EEXC=0.0;
EINH=-90.;//-75.0;
ZETA=0.5;
PHI=0.85;

// Type-II and IA=0
CM=0.012;// in nF
EK=-70.0;
ENA=55.0;
EH=-43.0;
ELK=-65.0;
GNA=1000.0;// in nS
GHT=150.0;
GLT=200.;//200.0;
GA=0.0;
GH=20.;//20.0;
GLK=2.;//2.0;

T=(1./(freq*1.));
tTOT=T*(Nspikes*1.+2.);
dt=T/Ndiscret;
ITDmax=1.5;
sigmaT1=.1*T;
sigmaT2=.1*T;
sigmaT3=.1*T;
sigmaT4=.1*T;
sigmaA1=.1*GEXC;
sigmaA2=.1*GEXC;
sigmaA3=.1*GINH1;
sigmaA4=.1*GINH2;

R1=0.;
R2=0.;
R3=0.;
R4=0.;

//************//
//Loop in ITDs//
//************//
o=0;
for(o=0;o<=Npoints;o++){
	//**************//
	//Loop in trials//
	//**************//
	for(p=0;p<N_steps_sta;p++){
		sta[p]=0.;
		sta_stdev[p]=0.;
		epsp_ave[p]=0.;
		epsp_stdev[p]=0.;
		uper[p]=-120.;
		epsp_ave_upper[p]=-120.;
		lower[p]=60.;
		epsp_ave_lower[p]=60.;
		iklt_ave[p]=0.;
		gklt_ave[p]=0.;
		epsp_gklt_ave[p]=0.;
		syn_ex1[p]=0.; 
		syn_ex2[p]=0.;
		syn_spike_ex1[p]=0.;
		syn_spike_ex2[p]=0.;
	}
	N_over_trials=0;
	Npeak=0;
	
	h=0;
	for(h=0;h<Ntrials;h++){
		//********************************//
		//Initianlization of the variables//
		//********************************//
		v[0]=-63.628399;
		v[1]=0.;
		v[2]=0.;
		v[3]=0.;
		v[4]=0.51221395;
		v[5]=0.66181272;
		v[6]=0.0077282726;
		v[7]=0.0011447769;
		v[8]=0.025057629;
		v[9]=0.44309759;
		v[10]=0.14586952;
		v[11]=0.;
		v[12]=0.;
		v[13]=0.;
		v[14]=0.;
		v[15]=0.;
		v[16]=0.;
		v[17]=0.;
		v[18]=0.;
		
		V0=v[0];
		w0=v[4];	
		z0=v[5];
		n0=v[6];
		p0=v[7];
		m0=v[8];
		h0=v[9];
		r0=v[10];
		
		N=0;
		P=0;
		v1=0.;
		v2=0.;
		v3=0.;
		for(p=0;p<N_steps_sta;p++){
			memory[p]=-63.628399;
			syn1[p]=0.;
			syn2[p]=0.;
			iklt_memory[p]=q3p03*GLT*(v[4]*v[4]*v[4]*v[4])*v[5]*(-63.628399-EK);
			gklt_memory[p]=q3p03*GLT*(v[4]*v[4]*v[4]*v[4])*v[5];
		}
		
		Nsmallepsp1=0;
		Nsmallepsp2=0;
		Nsmallepsp3=0;
		Nsmallepsp4=0;
		k=0;
		k1=0;
		k2=0;
		k3=0;
		k4=0;
		R1x=0.;
		R1y=0.;
		R2x=0.;
		R2y=0.;
		R3x=0.;
		R3y=0.;
		R4x=0.;
		R4y=0.;
		
		i=0;
		for(i=0;i<Nneurons_exc;i++){
			e1[i]=0.;
			e2[i]=0.;
			old1[i]=0;
			old2[i]=0;
		}
		
		i=0;
		for(i=0;i<Nneurons_inh;i++){
			e3[i]=0.;
			e4[i]=0.;
			old3[i]=0;
			old4[i]=0;
		}
		
		//*********************************************//
		//EPSPs and IPSPs time delay generation for ITD//		
		//*********************************************//				
		
		exc_ipsi_delay=o*ITDmax*T/(1.*Npoints);
		inh_ipsi_delay=o*ITDmax*T/(1.*Npoints)+0.2;
		inh_contra_delay=(ITDmax/2.)*T-0.6;
		exc_contra_delay=(ITDmax/2.)*T;
				
		//************//
		//Loop in time//
		//************//
		release_now=0.;
		t=0.;
		while(t<tTOT){
			
			//******************************************************************************//
			// lambda  //  nu  //  R=Vect_strength //  % of synapses active pre input-cycle // 
			// 0.75		  95       0.95					0.2									//	
			// 0.9		  24	   0.98					0.2									//
			// 0.55		  230      0.9					0.2									//
 			// 0.05		  800      0.8					0.2									//
			// -0.35	  1300     0.7					0.2									//	
			// -0.75	  1750     0.6					0.2									//	
			// -0.99	  2500     0.5					0.2									//	
			//
			//******************************************************************************//
		
			lambda1=0.975;
			lambda2=0.975;
			lambda3=0.975;
			lambda4=0.975;
			threshold1=0.9;//0.55;//0.75;//0.90;
			threshold2=-.99;//0.55;//0.75;//-0.99;
			threshold3=-.99;//-0.75;//0.9;
			threshold4=-.99;//-0.75;//0.9;//.75
			nu1=24;
			nu2=2500;
			nu3=2500;
			nu4=2500;
			
			prob_max_ipsi_exc=sin(2.*PI*(t - exc_ipsi_delay)*freq);
			prob_max_contra_exc=sin(2.*PI*(t - exc_contra_delay)*freq);
			prob_max_ipsi_inh=sin(2.*PI*(t - inh_ipsi_delay)*freq);
			prob_max_contra_inh=sin(2.*PI*(t - inh_contra_delay)*freq);
			if(t<exc_contra_delay){prob_max_ipsi_exc=-1.;}
			if(t<exc_ipsi_delay){prob_max_contra_exc=-1.;}
			if(t<inh_ipsi_delay){prob_max_ipsi_inh=-1.;}
			if(t<inh_contra_delay){prob_max_contra_inh=-1.;}
			
		
			//IPSILATERAL EPSPs
				i=0;
				for(i=0;i<Nneurons_exc;i++){
					release_now=(nu1*(2.*drand48())-1.);//24 for 0.90//95 for 0.75///230 for 0.55//800 for 0.05//1250 -0.35//1700 for -0.75// 2300 for -0.99 //2500 for -0.999 //  
					if((threshold1 < release_now )&&(release_now < prob_max_ipsi_exc)){
						if((e1[i])<t){
							v[12]=v[12]+GEXC1*exp(1.)/tauex1;
							old1[i]=e1[i];
							e1[i]=t;
							//fprintf(prnt4,"%lg\t%d\n",e1[i],i);
							k1=k1+1;
							Nsmallepsp1=k1;
							if(k1>1){
								R1x=R1x+cos(e1[i]*2.*PI*freq);
								R1y=R1y+sin(e1[i]*2.*PI*freq);
							}
						}
					}
				}
			
			//CONTRALATERAL EPSPs	
				i=0.;
				for(i=0;i<Nneurons_exc;i++){
					release_now=(nu2*(2.*drand48())-1.);
					if((threshold2 < release_now )&&(release_now < prob_max_contra_exc)){
						if((e2[i])<t){
							v[14]=v[14]+GEXC2*exp(1.)/tauex2;
							old2[i]=e2[i];
							e2[i]=t;
							//fprintf(prnt4,"%lg\t%d\n",e2[i],i);
							k2=k2+1;
							Nsmallepsp2=k2;
							if(k2>1){
								R2x=R2x+cos(e2[i]*2.*PI*freq);
								R2y=R2y+sin(e2[i]*2.*PI*freq);
							}
						}
					}
				}
			
			//IPSILATERAL IPSPs	
				i=0.;
				for(i=0;i<Nneurons_inh;i++){
					release_now=(nu3*(2.*drand48())-1.);
					if((threshold3 < release_now )&&(release_now < prob_max_ipsi_inh)){
						if((e3[i])<t){
							v[16]=v[16]+GINH1*exp(1.)/tauinh2;
							old3[i]=e3[i];
							e3[i]=t;
							//fprintf(prnt4,"%lg\t%d\n",e3[i],i);
							k3=k3+1;
							Nsmallepsp3=k3;
							if(k3>1){
								R3x=R3x+cos(e3[i]*2.*PI*freq);
								R3y=R3y+sin(e3[i]*2.*PI*freq);
							}
						}
					}
				}
						
			//CONTRALATERAL IPSPs	
			i=0.;
				for(i=0;i<Nneurons_inh;i++){
					release_now=(nu4*(2.*drand48())-1.);//80//1700//2200
					if((threshold4 < release_now )&&(release_now < prob_max_contra_inh)){
						if((e4[i])<t){
							v[18]=v[18]+GINH2*exp(1.)/tauinh2;
							old4[i]=e4[i];
							e4[i]=t;
							//fprintf(prnt4,"%lg\t%d\n",e4[i],i);
							k4=k4+1;
							Nsmallepsp4=k4;
							if(k4>1){
								R4x=R4x+cos(e4[i]*2.*PI*freq);
								R4y=R4y+sin(e4[i]*2.*PI*freq);
							}
						}
					}
				}	
				
			// Fast transient K+ current
			a_inf=1./(pow((1.0+exp(-(v[0]+31.0)/6.0)),0.25));
			b_inf=1./(pow((1.0+exp((v[0]+66.0)/7.0)),0.5));
			c_inf=b_inf;
			tau_a=q0p17*(100.0/(7.0*exp((v[0]+60.0)/14.0)+29.0*exp(-(v[0]+60.0)/24.0))+0.1);
			tau_b=q0p17*(1000.0/(14.0*exp((v[0]+60.0)/27.0)+29.0*exp(-(v[0]+60.0)/24.0))+1.0);
			tau_c=q0p17*(90.0/(1.0+exp(-(v[0]+66.0)/17.0))+10.0);
                                                                                            
			// Low threshold K+ current
			w_inf=1./(pow((1.0+exp(-(v[0]+48.0)/6.0)),0.25));
			z_inf=(1.0-ZETA)/(1.0+exp((v[0]+71.0)/10.0))+ZETA;
			tau_w=0.8*q0p17*(100.0/(6.0*exp((v[0]+60.0)/6.0)+16.0*exp(-(v[0]+60.0)/45.0))+1.5);
			tau_z=q0p17*(1000.0/(exp((v[0]+60.0)/20.0)+exp(-(v[0]+60.0)/8.0))+50.0);
                                                                                            
			// High threshold K+ current
			n_inf=1./(pow((1.0+exp(-(v[0]+15.0)/5.0)),0.5));
			p_inf=1.0/(1.0+exp(-(v[0]+23.0)/6.0));
			tau_n=q0p17*(100.0/(11.0*exp((v[0]+60.0)/24.0)+21.0*exp(-(v[0]+60.0)/23.0))+0.7);
			tau_p=q0p17*(100.0/(4.0*exp((v[0]+60.0)/32.0)+5.0*exp(-(v[0]+60.0)/22.0))+5.0);
                                                                                            
			// Fast Na+ current
			m_inf=1.0/(1.0+exp(-(v[0]+38.0)/7.0));
			h_inf=1.0/(1.0+exp((v[0]+65.0)/6.0));
			tau_m=q0p17*(10.0/(5.0*exp((v[0]+60.0)/18.0)+36.0*exp(-(v[0]+60.0)/25.0))+0.04);
			tau_h=q0p17*(100.0/(7.0*exp((v[0]+60.0)/11.0)+10.0*exp(-(v[0]+60.0)/25.0))+0.6);
                                                                                            
			// Hyperpolarization-activated cation current
			r_inf=1.0/(1.0+exp((v[0]+76.0)/7.0));
			tau_r=q0p17*(100000./(237.0*exp((v[0]+60.0)/12.0)+17.0*exp(-(v[0]+60.0)/14.0))+25.0);
                                                                                            
			// Current equations
			gA=0.0;//GA*(v[1]*v[1]*v[1]*v[1])*v[2]*v[3]*(v[0]-EK)
			gLT=1.3*q3p03*GLT*(v[4]*v[4]*v[4]*v[4])*v[5];//1.6*q3p03*GLT*(v[4]*v[4]*v[4]*v[4])*v[5];//1.6*q3p03*GLT*(w0*w0*w0*w0)*z0;//1.6*q3p03*GLT*(v[4]*v[4]*v[4]*v[4])*v[5];
			gHT=q3p03*GHT*(PHI*v[6]*v[6]+(1.0-PHI)*v[7]);
			gNA=q3p03*GNA*(v[8]*v[8]*v[8])*v[9];
			gH=1.3*q3p03*GH*v[10];//1.56*q3p03*GH*r0;//1.56*q3p03*GH*v[10];
			gLK=q3p03*GLK;
            
			IA=gA*(v[0]-EK);
			ILT=gLT*(v[0]-EK);
			IHT=gHT*(v[0]-EK);
			INA=gNA*(v[0]-ENA);
			IH=gH*(v[0]-EH);
			ILK=gLK*(v[0]-ELK);
                                                                                            
			cur_term=I0;
			
			syn1ex = v[11]*(v[0]-EEXC);
			syn2ex = v[13]*(v[0]-EEXC);
			syn1in = v[15]*(v[0]-EINH);
			syn2in = v[17]*(v[0]-EINH);
			syn_term= syn1ex + syn2ex + syn1in +syn2in;

			rk4(takens,v,19,t,dt);
			
			//******************************************************************************************************//
			//Loop to compute Spike triggered voltage, synaptic conductance, IKLT current, IKLT conductance averages//
			//******************************************************************************************************//
			
			
			for(p=0;p<(N_steps_sta-1);p++){
				memory[N_steps_sta - p -1] = memory[N_steps_sta - p - 2];
				syn1[N_steps_sta - p -1]=syn1[N_steps_sta - p - 2];
				syn2[N_steps_sta - p -1]=syn2[N_steps_sta - p - 2];
				iklt_memory[N_steps_sta - p -1] = iklt_memory[N_steps_sta - p - 2];
				gklt_memory[N_steps_sta - p -1] = gklt_memory[N_steps_sta - p - 2];
			}	
			memory[0]=v[0];
			iklt_memory[0]=ILT;
			gklt_memory[0]=q3p03*GLT*(v[4]*v[4]*v[4]*v[4])*v[5];
			syn1[0]=syn1ex;
			syn2[0]=syn2ex;
				
			v3=v2;
			v2=v1;
			v1=v[0];
			
			if(((v2-v3)>0.)&&((v1-v2)<0.)){
				P=P+1;
				p=0;
				for(p=0;p<N_steps_sta;p++){
					epsp_ave[p]=epsp_ave[p]+memory[N_steps_sta-p-1];
					epsp_stdev[p]=epsp_stdev[p]+memory[N_steps_sta-p-1]*memory[N_steps_sta-p-1];
					epsp_gklt_ave[p]=epsp_gklt_ave[p]+gklt_memory[N_steps_sta - p -1];
					syn_ex1[p]=syn_ex1[p] + syn1[N_steps_sta-p-1]; 
					syn_ex2[p]=syn_ex2[p] + syn2[N_steps_sta-p-1];
					if(memory[N_steps_sta-p-1]<epsp_ave_lower[p]){epsp_ave_lower[p]=memory[N_steps_sta-p-1];}
					if(epsp_ave_upper[p]<memory[N_steps_sta-p-1]){epsp_ave_upper[p]=memory[N_steps_sta-p-1];}
				}	
			}
			if(((v2-v3)>0.)&&((v1-v2)<0.)&&(v2>-25.)){
				N=N+1;
				p=0;
				for(p=0;p<N_steps_sta;p++){
					sta[p]=sta[p] + memory[N_steps_sta-p-1];
					sta_stdev[p]=sta_stdev[p] + memory[N_steps_sta-p-1]*memory[N_steps_sta-p-1];
					iklt_ave[p]=iklt_ave[p] + iklt_memory[N_steps_sta-p-1];
					gklt_ave[p]=gklt_ave[p] + gklt_memory[N_steps_sta-p-1];
					syn_spike_ex1[p]=syn_spike_ex1[p] + syn1[N_steps_sta-p-1]; 
					syn_spike_ex2[p]=syn_spike_ex2[p] + syn2[N_steps_sta-p-1];
					if(memory[N_steps_sta-p-1]<lower[p]){lower[p]=memory[N_steps_sta-p-1];}
					if(uper[p]<memory[N_steps_sta-p-1]){uper[p]=memory[N_steps_sta-p-1];}
				}
				
			}
			
			if(o==Npoints){
				fprintf(prnt2,"%lg\t%lg\t%lg\t%lg\t%lg\t%lg\n",t,v[0],sin(2.*PI*freq*(t)),syn_term,v[11],v[13]);
			}
			
			t=t+dt;
			
		}
		
		//***********************************//
		//END OF THE LOOP FOR THE RK INTEGRATOR
		//***********************************//
		
		N_over_trials=N_over_trials+N;
		Npeak=Npeak+P;
		
		R1 = R1+(sqrt(R1x*R1x+R1y*R1y)/(1.*(Nsmallepsp1-1)));
		R2 = R2+(sqrt(R2x*R2x+R2y*R2y)/(1.*(Nsmallepsp2-1)));
		R3 = R3+(sqrt(R3x*R3x+R3y*R3y)/(1.*(Nsmallepsp3-1)));
		R4 = R4+(sqrt(R4x*R4x+R4y*R4y)/(1.*(Nsmallepsp4-1)));

	}
	
	prnt1=fopen("Model1_STA_epsg_epsp_iklt_w-inh.dat","a");
	p=0;
	for(p=0;p<N_steps_sta;p++){
		fprintf(prnt1,"%lg\t%d\t%lg\t%lg\t%lg\t%lg\t%lg\t%lg\t%lg\t%lg\t%lg\t%lg\t%lg\t%lg\t%lg\t%lg\t%lg\n",exc_ipsi_delay-(ITDmax/2.)*T,p,1.*sta[p]/(1.*N_over_trials),lower[p],uper[p],epsp_ave[p]/(1.*Npeak),epsp_ave_lower[p],epsp_ave_upper[p],sqrt((epsp_stdev[p]/(1.*Npeak))-(epsp_ave[p]/(1.*Npeak))*(epsp_ave[p]/(1.*Npeak))),sqrt((sta_stdev[p]/(1.*N_over_trials))-(sta[p]/(1.*N_over_trials))*(sta[p]/(1.*N_over_trials))),iklt_ave[p]/(1.*N_over_trials),gklt_ave[p]/(1.*N_over_trials),epsp_gklt_ave[p]/(1.*Npeak),syn_ex1[p]/(1.*Npeak),syn_ex2[p]/(1.*Npeak),syn_spike_ex1[p]/(1.*N_over_trials),syn_spike_ex2[p]/(1.*N_over_trials));
	}
	
	fprintf(prnt1,"\n");	
	fclose(prnt1);
	
	prnt3=fopen("Model1_ITD-curve_w-inh.dat","a");
	fprintf(prnt3,"%lg\t%lg\n",exc_ipsi_delay-(ITDmax/2.)*T,(1.*N_over_trials)/(1.*Ntrials));
	fclose(prnt3);
	
}	
	    
    R1 = sqrt( (R1x/(1.*Nsmallepsp1-1.))*(R1x/(1.*Nsmallepsp1-1.)) + (R1y/(1.*Nsmallepsp1-1.))*(R1y/(1.*Nsmallepsp1-1.)) );
    R2 = sqrt( (R2x/(1.*Nsmallepsp2-1.))*(R2x/(1.*Nsmallepsp2-1.)) + (R2y/(1.*Nsmallepsp2-1.))*(R2y/(1.*Nsmallepsp2-1.)) );
    R3 = sqrt( (R3x/(1.*Nsmallepsp3-1.))*(R3x/(1.*Nsmallepsp3-1.)) + (R3y/(1.*Nsmallepsp3-1.))*(R3y/(1.*Nsmallepsp3-1.)) );		
    R4 = sqrt( (R4x/(1.*Nsmallepsp4-1.))*(R4x/(1.*Nsmallepsp4-1.)) + (R4y/(1.*Nsmallepsp4-1.))*(R4y/(1.*Nsmallepsp4-1.)) );
	 
	printf("Exc Ipsi Vect-Stren= ");
	printf("%lg\t%lg\n",R1,(1.*Nsmallepsp1)/(1.*Nneurons_exc)/Nspikes);
	printf("Exc Contra Vect-Stren= ");
    printf("%lg\t%lg\n",R2,(1.*Nsmallepsp2)/(1.*Nneurons_exc)/Nspikes);
	printf("Inh Ipsi Vect-Stren= ");
	;printf("%lg\t%lg\n",R3,(1.*Nsmallepsp3)/(1.*Nneurons_inh)/Nspikes);
	printf("Inh Contra Vect-Stren= ");
	printf("%lg\t%lg\n",R4,(1.*Nsmallepsp4)/(1.*Nneurons_inh)/Nspikes);		

fclose(prnt1);
fclose(prnt2);
fclose(prnt3);

Gtotal=(q3p03*GLT*(w0*w0*w0*w0)*z0+q3p03*GHT*(PHI*n0*n0+(1.0-PHI)*p0)+q3p03*GNA*(m0*m0*m0)*h0+q3p03*GH*r0+q3p03*GLK); 

printf("membrane time constant= ");
printf("%lg",(CM/Gtotal)*1000. );
printf("msec \n");
printf("membrane resistance= ");
printf("%lg",(1./Gtotal));
printf("Ohm/cm^2 \n");

exit(0);

}

Loading data, please wait...