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 32
#define PI 3.1415926535897932384626433832795028841971 
#define Ntrials 5
#define Npoints 20
#define Nspikes 25
#define TRUE 1
#define FALSE 0
#define Ndiscret 1000
#define N_steps_sta 1000

double t_q,R_r,R_st,t_r,R_st,t_st,R_ss,phi,phi1,phi2,phi01,phi02,freq,nu,nu1,nu2,tau_d,t_k; 
double Rss,RA,C0,C1,S0,S1;
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;

/* Neuron model with synaptic inputs equations */

void takens(int n, double v[], double dv[], double t)
{
 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 - 1.*v[11]/(tauex1*tauex1);
dv[13]=v[14];
dv[14]= - 2.*v[14]/tauex2 - 1.*v[13]/(tauex2*tauex2);
dv[15]=v[16];
dv[16]= - 2.*v[16]/tauinh1 - 1.*v[15]/(tauinh1*tauinh1);
dv[17]=v[18];
dv[18]= - 2.*v[18]/tauinh2 - 1.*v[17]/(tauinh2*tauinh2);

return;
}


/* This routine computes the nth stage of refinement of an extended trapezoidal rule. func is input as a pointer to the function to be 
 integrated between limits a and b, also input. When called with	n=1, the routine returns the crudest estimate of integral from a 
 to b of f(x)dx. Subsequent calls with n=2,3,...(in that sequential order) will improve the accuracy by adding 2n-2 additional interior 
 points. */

#define FUNC(x) ((*func)(x))

double trapzd(double (*func)(double), double a, double b, int n){		

	double x,tnm,sum,del;
	static double s;
	int it,j;
	if (n == 1) {
		return (s=0.5*(b-a)*(FUNC(a)+FUNC(b)));
	} else {
	for (it=1,j=1;j<n-1;j++) it <<= 1;
		tnm=it;
		del=(b-a)/tnm; 						//This is the spacing of the points to be added.
		x=a+0.5*del;
		for (sum=0.0,j=1;j<=it;j++,x+=del) sum += FUNC(x);
		s=0.5*(s+(b-a)*sum/tnm); 				//This replaces s by its refined value.
		return s;
	}
}


/* This routine returns the integral of the function func from a to b. The parameters EPS can be set to the desired fractional accuracy and 
 JMAX so that 2 to the power JMAX-1 is the maximum allowed number of steps. Integration is performed by the trapezoidal rule. */ 

#define EPS 1.0e-3
#define JMAX 20
double qtrap(double (*func)(double), double a, double b){ 
double trapzd(double (*func)(double), double a, double b, int n);
int j;
double s,olds=0.0; 					//Initial value of olds is arbitrary.

	for (j=1;j<=JMAX;j++) {
		s=trapzd(func,a,b,j);
		if (j > 5) 				//Avoid spurious early convergence.
		if (fabs(s-olds) < EPS*fabs(olds)||(s == 0.0 && olds == 0.0)) return s;
		olds=s;
	}
	return 0.0; 					//Never get here.
}

/*#############################################################################################################
Definition of the spontaneuos figing rate functions to be evaluated by a time dependent poission process scheme. 
##############################################################################################################*/ 

double lambdaIpsi(double t){
	return(( R_ss*(exp(phi1*sin(2.*PI*freq*(t)))/exp(phi1))+phi01)*nu1*(1.-0.55*exp(-(t-t_k)/(0.8))-0.45*exp(-(t-t_k)/(25.)))); //0.0008 & 0.09 //Vs=0.88: 0.0025 & 0.13
}

double lambdaContra(double t){
	return(( R_ss*(exp(phi2*sin(2.*PI*freq*(t)))/exp(phi2))+phi02)*nu2*(1.-0.55*exp(-(t-t_k)/(0.8))-0.45*exp(-(t-t_k)/(25.))));//*(1.-exp(nu*(t-t-1))));//*C0*exp(-(t)/S0)); //0.0008 & 0.09 //Vs=0.88: 0.0025 & 0.13
}

/*Other functions:
double lambda(double t){
	return( (1.-exp(-(t+t_k)/t_q))*( R_r*exp(-(t+t_k)/t_r) + R_st*exp(-(t+t_k)/t_st) + R_ss*(exp(phi*sin(2.*PI*freq*(t+t_k)))/exp(phi)))*(1.-exp(-nu*(t))) ); //0.0008 & 0.09 //Vs=0.88: 0.0025 & 0.13
}
double lambda1(double t){
	return( R_ss*(exp(sin(2.*PI*freq*(t+t_k)))/(exp(phi)))*(1.-exp(-nu*(t))));
}
*/


/*########################################
 Local variables definition and main MSO model code
######################################## */

int main(){
FILE *prnt0, *prnt1, *prnt2, *prnt3, *prnt4,*prnt5;
double v[19], e1[Nneurons], e2[Nneurons],e3[Nneurons],old1[Nneurons],old2[Nneurons],old3[Nneurons];
double I0, IEXT,q0p17,q3p03;
double startex,EEXC,startinh,EINH,ZETA,PHI;
double EK,ENA,EH,ELK,GNA,GHT,GLT,GA,GH,GLK,VREST,RREST,TAUM,VTH,SM,GE22,GE38;
double t,t_int,dt,tTOT,T,ITDmax;
double sigma1,sigma2,sigma3,d1,d2,d3,v1,v2,v3,integral;
double u,u1,u2,R,Phi,nu0;
int i,j,k,l,o,h,aa,pp;
int N,decision,Npeak,N_over_trials;
int Nsmallepsp1,Nsmallepsp2,Nsmallepsp3,k1,k2,k3;
double tis,R1x,R1y,R2x,R2y,R3x,R3y,R1,R2,R3,y2;
double x[Nneurons],y[Nneurons],yy[Nneurons];
double Gtotal, gA,gHT,gLT,gNA,gH,gLK;
double V0,w0,z0,n0,p0,m0,h0,r0;	
double max,twenty,eighty,t1,t2,t3,t4,peak,v0;	
double aver_epsp[Ndiscret],aver_epsg_ips[Ndiscret],aver_epsg_con[Ndiscret],aver_ipsg[Ndiscret],sta[N_steps_sta],memory[N_steps_sta];

//#######################	
//Output Files definition
//#######################
	
prnt1=fopen("Model2_syn-events-ipsi.dat","wt");
prnt2=fopen("Model2_syn-events-contra.dat","wt");	
prnt3=fopen("Model2_train_EPSGs-EPSPs.dat","wt");
prnt4=fopen("Model2_ITD-curve.dat","w");
prnt5=fopen("Model2_STA_epsg_epsp_iklt.dat","wt");

	/*Paramter values for lambda functions*/
	
	freq=.5;//in KHz
	phi=7.70 - 3.69*freq + 0.40*freq*freq;//3.79 - 1.08*freq + 0.07*freq*freq;
	t_q=0.2;
	t_r=3.;
	t_st=10.;
	R_r=600.;
	R_st=200.;
	R_ss=2296 - 977*freq + 105*freq*freq;
	
	nu=(10./8.)*(R_r+R_st+2.*R_ss);
	
	nu1=0.006;//0.0065;//0.0018;
	nu2	=0.0009;//0.0042;//0.0005;//0.001;
	phi1=32.7525;//7.8606;//phi*1.45;//phi1=phi*10.,nu1=0.01875,VS=0.879479,%=1.00875,Slope=4.82404,HW=0.586,GEXC=1.15//phi1=phi*10.,nu1=0.003,VS=0.889077,%=0.34125,Slope=4.92063,HW=0.602
	phi2=1.5483;//1.7865;//1.0719;//phi*.25;//phi2=phi*.1,nu2=0.00105,VS=0.236316,%=1.03,Slope=1.72839,HW=0.32//phi2=phi*.25,nu2=0.00185,VS=0.482553,%=0.97125,Slope=2.09915,HW=0.422//phi2=phi*.2,nu2=0.00035;VS=0.483799,%=0.35125,Slope=1.64128,HW=0.436
	phi01=0.0;//-0.05;
	phi02=0.;//-0.1;//0.0;	
	
//	32.7525 1.7865  0       0       0.005   0.0042  0.196   -0.24
	
	
//#######################
//Experimental setting for MSO neuron cell    
//#######################
I0=0.0;
IEXT=0.0;
q0p17=0.17;//1.0;
q3p03=3.03;//1.0;
	

//#######################	
//parameters for syn_term:
//#######################	
	GEXC1=3.;//3.;//2.5;//.9;//.5;//.50;//1.35;
	GEXC2=4.5;//4.5;//4.8;//.5;//.50;//1.35;
GINH1=0.0;//.70;
GINH2=0.0;
syn1ex=0.;
syn2ex=0.;
syn1in=0.;
syn2in=0.;
tauex1=0.10;
tauex2=0.10;
tauinh1=.10;
tauinh2=.10;
EEXC=0.0;
EINH=-90.;//-75.;
ZETA=0.5;
PHI=0.85;
                                                                                            
// Type-II and IA=0
CM=0.012;
EK=-70.0;
ENA=55.0;
EH=-43.0;
ELK=-65.0;
GNA=1000.0;
GHT=150.0;
GLT=200.0;
GA=0.0;
GH=20.0;
GLK=2.0;

//#######################	
//Simulations time values.
//#######################	
T=(1./(freq*1.));
tau_d=.75;
tTOT=T*(Nspikes*1.+2.);
dt=T/Ndiscret;
ITDmax=.6;	
srand48(time(0));	

    o=0;
    for(o=0;o<=Npoints;o++){
		l=0;
		N_over_trials=0;
		for(l=0;l<Ntrials;l++){
			
			//#####################################
			//Inizialization of dynamical 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;
			v1=0.;
			v2=0.;
			v3=0.;
		
			d1=o*ITDmax*T/(1.*Npoints);
			d2=(ITDmax/2.)*T;
			d3=0.0;		

			Nsmallepsp1=0;
			Nsmallepsp2=0;
			Nsmallepsp3=0;
			k=0;
			k1=0;
			k2=0;
			k3=0;
			R1x=0.;
			R1y=0.;
			R2x=0.;
			R2y=0.;
			R3x=0.;
			R3y=0.;
		
			i=0;
			for(i=0;i<Nneurons;i++){
        		e1[i]=drand48()*(T/5.);
        		e2[i]=drand48()*(T/5.);
        		e3[i]=drand48()*(T/5.);
        		old1[i]=0.;
        		old2[i]=0.;
        		old3[i]=0.;
			}
			
			for(aa=0;aa<Ndiscret;aa++){
				aver_epsp[aa]=v[0];
				aver_epsg_ips[aa]=v[11];
				aver_epsg_con[aa]=v[13];
				aver_ipsg[aa]=v[17];
			}
			
			t=0.;
			aa=0;

			//%%%%%%//Main simulation time loop//%%%%%%%%%
			
			while(t<tTOT){
				Npeak=floor(t/T);
				
				if(fabs(fmod(t,T)-2.)<dt)aa=0;
				
				if(o==Npoints){
					aver_epsp[aa]=aver_epsp[aa]+v[0];
					aver_epsg_ips[aa]=aver_epsg_ips[aa]+v[11];
					aver_epsg_con[aa]=aver_epsg_con[aa]+v[13];
					aver_ipsg[aa]=aver_ipsg[aa]+v[17];
				}
			
				//Decision to make an event on the ipsilateral input
				i=0;
				for(i=0;i<Nneurons;i++){
					decision=TRUE;
					if((fabs(e1[i]-t+d1)<0.5*dt)){
						u=drand48();
						t_k=e1[i];
						t_int=tau_d;
						decision=FALSE;
						while(decision == FALSE){
							if(t_int>tTOT)decision=TRUE;
							
							integral=qtrap(&lambdaIpsi,t_k,t_k+t_int);
							
							if(((log(u)+integral)>=0.)){
								old1[i]=e1[i];
								v[12]=v[12]+GEXC1*exp(1.)/tauex1;
								e1[i]=t_int+old1[i];
								fprintf(prnt1,"%lg\t%d\t%lg\t%lg\n",t,i,e1[i],lambdaIpsi);
								decision=TRUE;
							}
							
							t_int=t_int+10.*dt;
							
						}
						
						k1=k1+1;
						Nsmallepsp1=k1;
						
						if(k1>1){
							R1x=R1x+cos(e1[i]*2.*PI*freq);
							R1y=R1y+sin(e1[i]*2.*PI*freq);
						}
						
					}
				}
				
				//Decision to make an event on the contralateral input
				j=0;
				for(j=0;j<Nneurons;j++){
					decision=TRUE;
					if((fabs(e2[j]-t+d2)<0.5*dt)){
						u=drand48();
						t_k=e2[j];
						t_int=tau_d;
						decision=FALSE;
						while(decision == FALSE){
							if(t_int>tTOT)decision=TRUE;
							
							integral=qtrap(&lambdaContra,t_k,t_k+t_int);
							
							if((log(u)+integral)>=0. ){
								old2[j]=e2[j];
								v[14]=v[14]+GEXC2*exp(1.)/tauex2;
								e2[j]=t_int+old2[j];
								fprintf(prnt2,"%lg\t%d\t%lg\t%lg\n",t,j,e2[j],lambdaContra);
								decision=TRUE;
							}
							
							t_int=t_int+10.*dt;
							
						}
						k2=k2+1;
						Nsmallepsp2=k2;
						
						if(k2>1){
							R2x=R2x+cos(e2[j]*freq*2.*PI);
							R2y=R2y+sin(e2[j]*freq*2.*PI);
						}
					}
				}


					//#######################
					//MSO Point Neuron Model	
					//#######################
					// 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=q0p17*(90.0/(6.0*exp((v[0]+60.0)/6.0)+16.0*exp(-(v[0]+60.0)/45.0))+0.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
					IA=0.0;//GA*(v[1]*v[1]*v[1]*v[1])*v[2]*v[3]*(v[0]-EK)
					ILT=1.3*q3p03*GLT*(v[4]*v[4]*v[4]*v[4])*v[5]*(v[0]-EK);//q3p03*GLT*(w0*w0*w0*w0)*z0*(v[0]-EK);//q3p03*GLT*(v[4]*v[4]*v[4]*v[4])*v[5]*(v[0]-EK);
					IHT=q3p03*GHT*(PHI*v[6]*v[6]+(1.0-PHI)*v[7])*(v[0]-EK);
					INA=q3p03*GNA*(v[8]*v[8]*v[8])*v[9]*(v[0]-ENA);
					IH=1.29*q3p03*GH*v[10]*(v[0]-EH);
					ILK=q3p03*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);
			
				//####################//Creating the spike triggered average//################
				pp=0;
				for(pp=0;pp<N_steps_sta;pp++){
					memory[N_steps_sta-pp]=memory[N_steps_sta - pp - 1];
				}	
				memory[1]=v[0];
				
				v3=v2;
				v2=v1;
				v1=v[0];
				if(((v2-v3)>0.)&&((v1-v2)<0.)&&(v2>-20.)){
					N=N+1;
					printf("%d\n",N);
					pp=0;
					for(pp=0;pp<N_steps_sta;pp++){
						sta[pp]=sta[pp]+memory[N_steps_sta-pp];
					}
				}
				
				
				if(o==Npoints){
					fprintf(prnt3,"%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;
				aa=aa+1;

		}
		fprintf(prnt2,"\n");
		
		//####################//Creating the Vector strength for each input//################
			
		N_over_trials=N_over_trials+N;	
			
		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)));

	
	}
	
		prnt3=fopen("Model2_ITD-curve.dat","a");
		fprintf(prnt4,"%lg\t%lg\n",d1-(ITDmax/2.)*T,((1.*N_over_trials)/(1.*Ntrials))/Nspikes);
		fclose(prnt4);	
		}
		
		//####################//Creating the average slope and halfwidth//################
	
		//REMEMBER TO PUT A "-" SIGN IN FRONT OF THE VECTOR VALUE IN CASE OF IPSP
		
		v0=(aver_epsp[1]/((1.*Nspikes+3.)));//-(aver_epsp[1]/((1.*Nspikes+3.)));
		v3=v2=v1=v0;
		max=0.;
		for(aa=1;aa<(Ndiscret-1);aa++){
			fprintf(prnt5,"%lg\t%lg\t%lg\t%lg\t%lg\n",aa*dt,aver_epsp[aa]/(1.*Nspikes+3.),aver_epsg_ips[aa]/(1.*Nspikes+3.),aver_epsg_con[aa]/(1.*Nspikes+3.),aver_ipsg[aa]/(1.*Nspikes+3.));
			v3=v2;
			v2=v1;
			v1=(aver_epsp[aa]/((1.*Nspikes+3.)));//-(aver_epsp[aa]/((1.*Nspikes+3.)));
			if((v3<v2)&&(v1<v2)&&(max<(v2-v0))){
				peak=aa*dt;
				max=v2-v0;
			}
		}
		
		for(aa=1;aa<(Ndiscret-1);aa++){
			v3=v2;
			v2=v1;
			v1=(aver_epsp[aa]/((1.*Nspikes+3.)));//-(aver_epsp[aa]/((1.*Nspikes+3.)));
			if(((v3-v0)<(max/2.))&&((max/2.)<(v2-v0))&&(aa*dt<(peak))){t1=aa*dt;}
			if(((v2-v0)<(max/2.))&&((max/2.)<(v3-v0))&&(aa*dt<(peak))){t2=aa*dt;}
			if(((v3-v0)<(max*0.2))&&((max*0.2)<(v2-v0))&&(aa*dt<(peak))){t3=aa*dt;twenty=v2-v0;}
			if(((v3-v0)<(max*0.8))&&((max*0.8)<(v2-v0))&&(aa*dt<(peak))){t4=aa*dt;eighty=v2-v0;}	
		}		
    
    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.)) );	
	
	//####################//Saving all the results on terminal//################
    
	printf("		Phi	R_ss	VS	Prcnt.	\n");
	printf("ipsi_exc	");
    printf("%lg\t%lg\t%lg\t%lg\n",phi1,R_ss*nu1,R1,(1.*Nsmallepsp1)/(1.*Nneurons)/Nspikes);
	printf("cont_exc	");
    printf("%lg\t%lg\t%lg\t%lg\n",phi2,R_ss*nu2,R2,(1.*Nsmallepsp2)/(1.*Nneurons)/Nspikes);
	printf("cont_inh	");
    printf("%lg\t%lg\t%lg\t%lg\n",phi,R_ss,R3,(1.*Nsmallepsp3)/(1.*Nneurons)/Nspikes);		
	
	printf("\n");
	
	fclose(prnt1);
	fclose(prnt2);
	fclose(prnt3);
	fclose(prnt4);
	fclose(prnt5);
	
	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("\n");
	printf("Slope	Norm-Slope	Halfwidth	\n");
	printf("%lg\t%lg\t%lg\n",(eighty-twenty)/(t4-t3),((eighty-twenty)/(t4-t3))/max,t1-t2);


exit(0);

}




Loading data, please wait...