Okamoto_etal_program.c

/*LINEAR ACCUMULATOR */
#include "math.h"
#include <stdio.h>
#include <stdlib.h>
#define L 505 /* cell no. */

/* parameters */
int gg;
int TRIAL;
double dt,time;
int Np,Ni; /* no. of neurons */
int CONNECTIONpp[L][L];
double connection_rate;
double Gpp,Gpi,Gip;

/* Membrane properties of a pyramidal cell */
double V_resting;
double V_theta;
double V_reset;
double G_L;
double refractory;
double Cp;

/* Membrane properties of an inter-neuron */
double Vi_resting;
double Vi_theta;
double Vi_reset;
double G_i_L;
double refractory_i;
double Ci;

/* Membrane potentials and related variables */
double V1[L],V2[L];
double Vi1[L],Vi2[L];
int AP[L];
int APi[L];
double t_last_AP[L];	/* lapse of time after the last AP */
double t_last_APi[L];	/* lapse of time after the last APi */

/* AMPA & GABA */
double E_AMPA,E_GABA;
double tau_AMPA,tau_GABA;
// Sparse connection
double p_AMPA[L],p_GABA[L];
double p_AMPA_jump,p_GABA_jump;

/* Population activity */
double FEFp,FEFi,FREQp,FREQi;

/* Recurrent */
double Ppp,Ppi,Pip;
// Sparse connection
double Ppp_RECURRENT[L];

/* ADP current */
double P_ADP[L],Ca[L],tau_Ca,K_Ca,n_Ca,del_Ca,E_ADP,G_ADP;

/* External current */
double I_base,I_base_i;
double I_dep,I_dep_value,time_dep_on,time_dep_off,I_hyp,I_hyp_value,I_hyp_on,I_hyp_off;

/* Noise */
double G_excitation[L],G_inhibition[L];
double G_excitation_i[L], G_inhibition_i[L];
double F_E,g_E,F_I,g_I;
double F_i_E,g_i_E,F_i_I,g_i_I;


/* Files */
FILE *fp_FEFp,*fp_FEFi;


/* RANDOM */
double RAND()
{
	double b;
	b=(double)rand()/(double)RAND_MAX;
	return(b);
}


/* initialization1 */
initialization1()
{
	char s[80];
	double RAND();
	double pp;
	int n,nn,m,mm,t;

	dt=0.00001;
	TRIAL=1;
	
	/* No. of pyramidal/inhibitory cells */
	Np=100;
	Ni=(int)(0.2*(double)Np);

	/* Membrane propertues of a pyramidal cell */
	V_resting=0.0-70.0;
	V_theta=0.0-52.0;
	V_reset=0.0-62.0;
	Cp=1.0;// scaled to 1.0
	G_L=(1.0/0.02)*Cp;
	refractory=0.002*2.0; 


	/* Membrane propertues of an inter-neuron */
	Vi_resting=0.0-65.0;
	Vi_theta=0.0-52.0;
	Vi_reset=0.0-60.0;
	Ci=1.0;// scaled to 1.0
	G_i_L=(1.0/0.01)*Ci;
	refractory_i=0.001*2.0;


	/* AMPA & GABA */
	E_AMPA=0.0;
	E_GABA=0.0-80.0;
	tau_AMPA=0.005;
	tau_GABA=0.005;
	p_AMPA_jump=0.50;
	p_GABA_jump=0.50;

	/* ADP */
	tau_Ca=0.07;
	K_Ca=1.0;
	n_Ca=4.0;
	del_Ca=1.3;
	E_ADP=0.0-35.0;
	G_ADP=12.0*Cp;

	/* Synaptic connection (conductance) */
	Gpp=20.0/(double)Np*Cp;	/* p->p */
	Gip=4.0/(double)Ni*Ci;		/* p->i */
	Gpi=20.0/(double)Np*Cp;	/* i->p */
	connection_rate=0.1;
		
	/* Noise in a pyramidal cell */
	F_E=(dt/1.0)*(1800.0);
	F_I=F_E;
	g_E=1.4*Cp;
	g_I=g_E*Cp;

	/* Noise in an interneuron */
	F_i_E=(dt/1.0)*(1800.0);
	F_i_I=F_i_E;
	g_i_E=1.4*Ci;
	g_i_I=g_i_E*Ci;

	/* Depolarization/hyperpolarization during the delay */
	// Depolarization
	I_dep_value=50.0*Cp;
	time_dep_on=0.0;
	time_dep_off=100.0;
	// Hyperpolarization
	I_hyp_value=(0.0-0.0*200.0)*Ci;
	I_hyp_on=4.0;
	I_hyp_off=4.2;

	/* Base current */
	I_base=400.0*Cp;
	I_base_i=1000.0*Ci;

	/* Connection matrix */
	for(n=1;n<=Np;++n)
	{
		for(nn=1;nn<=Np;++nn)
		{
			CONNECTIONpp[n][nn]=0;
			pp=RAND();
			if(pp<connection_rate)
			{
				CONNECTIONpp[n][nn]=1;
			}
		}
	}
	for(n=1;n<=Np;++n)
	{
		CONNECTIONpp[n][n]=0;
	}

}


/* initialization2 */
initialization2()
{
	char s[80];
	double RAND();
	int n;
	double pp;

	for(n=1;n<=Np;++n)
	{
		V2[n]=V_resting;
		t_last_AP[n]=10000.0;
		p_AMPA[n]=0.0;
		Ca[n]=0.0;
		G_excitation[n]=0.0;
		G_inhibition[n]=0.0;
	}

	for(n=1;n<=Ni;++n)
	{
		Vi2[n]=Vi_resting;
		t_last_APi[n]=10000.0;
		p_GABA[n]=0.0;
		G_excitation_i[n]=0.0;
		G_inhibition_i[n]=0.0;
	}

	FREQp=0.0;
	FREQi=0.0;

}

rewrite()
{
	int n,nn,m,mm;
	double RAND();

	for(n=1;n<=Np;++n)
	{
		V1[n]=V2[n];
		t_last_AP[n]=t_last_AP[n]+dt;
	}

	for(n=1;n<=Ni;++n)
	{
		Vi1[n]=Vi2[n];
		t_last_APi[n]=t_last_APi[n]+dt;
	}

	Ppp=0.0;

	for(n=1;n<=Np;++n)
	{
/*
		Ppp_RECURRENT[n]=0.0;
		for(nn=1;nn<=Np;++nn)
		{
			Ppp_RECURRENT[n]=Ppp_RECURRENT[n]
			+(double)CONNECTIONpp[n][nn]*p_AMPA[nn]*(1.0/connection_rate);
		}
*/
		Ppp=Ppp+p_AMPA[n];
		P_ADP[n]=pow(Ca[n],n_Ca)/(pow(Ca[n],n_Ca)+pow(K_Ca,n_Ca));
			
	}

	Pip=0.0;
	for(n=1;n<=Np;++n)
	{
		Pip=Pip+p_AMPA[n];
	}

	Ppi=0.0;
	for(n=1;n<=Ni;++n)
	{
		Ppi=Ppi+p_GABA[n];
	}

}


/* update */
update()
{
	int n,m,nn,mm;
	char s[80];
	double gauss();
	double RAND();
	double pp,pp_e,pp_i;

	/* Population activity */
	FREQp=FREQp*(1.0-dt/tau_AMPA);
	FREQi=FREQi*(1.0-dt/tau_GABA);
	FEFp=FREQp/(double)Np;
	FEFi=FREQi/(double)Ni;


	/* Excitatory cells */
	for(n=1;n<=Np;++n)
	{
		p_AMPA[n]=p_AMPA[n]*(1.0-dt/tau_AMPA);
		Ca[n]=Ca[n]*(1.0-dt/tau_Ca);
	
		pp_e=RAND();
		if(pp_e<F_E)
		{
			G_excitation[n]=G_excitation[n]+g_E;
		}
		G_excitation[n]=G_excitation[n]*(1.0-dt/tau_AMPA);
		pp_i=RAND();
		if(pp_i<F_I)
		{
			G_inhibition[n]=G_inhibition[n]+g_I;
		}
		G_inhibition[n]=G_inhibition[n]*(1.0-dt/tau_GABA);

		if(t_last_AP[n]>refractory)
		{
			V2[n]=V1[n]+dt*
			(
			(-1.0)*G_L*(V1[n]-V_resting)/Cp
			/*
			+Gpp*Ppp_RECURRENT[n]*(0.0-V1[n])/Cp
			*/
			+Gpp*(Ppp-p_AMPA[n])*(E_AMPA-V1[n])/Cp
			+Gpi*Ppi*(E_GABA-V1[n])/Cp
			+G_ADP*(E_ADP-V1[n])*P_ADP[n]/Cp
			+I_base/Cp
			+G_excitation[n]*(E_AMPA-V1[n])/Cp
			+G_inhibition[n]*(E_GABA-V1[n])/Cp
			+I_dep/Cp
			+I_hyp/Cp
			);
		}
		else
		{
			V2[n]=V1[n];
		}

		AP[n]=0;
		if(V2[n]>V_theta)
		{
			AP[n]=1;
			p_AMPA[n]=p_AMPA[n]+p_AMPA_jump*(1.0-p_AMPA[n]);
			t_last_AP[n]=0.0;
			/**/	
			FREQp=FREQp+1.0;
			V2[n]=V_reset;
			Ca[n]=Ca[n]+del_Ca;
		}
	}


	/* Inhibitory cells */	
	for(n=1;n<=Ni;++n)
	{
		p_GABA[n]=p_GABA[n]*(1.0-dt/tau_GABA);
		pp_e=RAND();
		if(pp_e<F_i_E)
		{
			G_excitation_i[n]=G_excitation_i[n]+g_i_E;
		}
		G_excitation_i[n]=G_excitation_i[n]*(1.0-dt/tau_AMPA);
		pp_i=RAND();
		if(pp_i<F_i_I)
		{
			G_inhibition_i[n]=G_inhibition_i[n]+g_i_I;
		}
		G_inhibition_i[n]=G_inhibition_i[n]*(1.0-dt/tau_GABA);

		if(t_last_APi[n]>refractory_i)
		{
			Vi2[n]=Vi1[n]+dt*
			(
			(-1.0)*G_i_L*(Vi1[n]-Vi_resting)/Ci
			+Gip*Pip*(E_AMPA-Vi1[n])/Ci
			+I_base_i/Ci
			+G_excitation_i[n]*(E_AMPA-Vi1[n])/Ci
			+G_inhibition_i[n]*(E_GABA-Vi1[n])/Ci
			);
		}
		else
		{
			Vi2[n]=Vi1[n];
		}

		APi[n]=0;
		if(Vi2[n]>Vi_theta)
		{
			APi[n]=1;
			p_GABA[n]=p_GABA[n]+p_GABA_jump*(1.0-p_GABA[n]);
			FREQi=FREQi+1.0;
			Vi2[n]=Vi_reset;
			t_last_APi[n]=0.0;
		}
	}

}

/* insert blanck line */
blanck_line()
{
	fprintf(fp_FEFp,"\n");
	fprintf(fp_FEFi,"\n");
}


/* file_print1 */
file_fprint1()
{
	fprintf(fp_FEFp,"%lf %lf\n",time,FEFp);
	fprintf(fp_FEFi,"%lf %lf\n",time,FEFi);
}


I_on_off()
{
	int n,nn,m,mm;
	double RAND();

		I_dep=0.0;
		if((time>time_dep_on)&&(time<time_dep_off))
		{
			I_dep=I_dep_value;
		}

		I_hyp=0.0;
		if((time>I_hyp_on)&&(time<I_hyp_off))
		{
			I_hyp=I_hyp_value;
		}
}



/* main program */ main() 
{
	char s[80];
	int R_SEED;

	/* file open */
	if((fp_FEFp=fopen("f_FEFp.out","a"))==NULL)
	{printf("File f_FEFp.out not exist\n");}
	if((fp_FEFi=fopen("f_FEFi.out","a"))==NULL)
	{printf("File f_FEFi.out not exist\n");}

	/* set a seed for random variables */
	printf("R_SEED:");
	sscanf(gets(s),"%d",&R_SEED);
	printf("\n");
	srand(R_SEED);

	initialization1();

	for(gg=1;gg<=TRIAL;++gg)
	{
		initialization2();
		time=0.0-2.0;

		blanck_line();
	//	system("del f_scatter.out");

		while(time<3.0)
		{
			time=time+dt;
			rewrite();
			I_on_off();
			update();
			if(((int)(time/dt)%50)==0)
			{
				file_fprint1();
			}
		}


		if((gg%1)==0)
		{
			printf("trial=%d\n",gg);
		}


	}

	/* file close */
	fclose(fp_FEFp);
	fclose(fp_FEFi);

	printf("Program terminated\n");
}