//Rulkov, N., Timofeev, I. and Bazhenov, M. "Oscillations in large-scale cortical //networks: map-based model" (2004) Journal of Computational Neuroscience 17, 203-224, 2004. #include #include #include #include //---------------Parameters initialization--------------------------------------------------- //--------------- Network Geometry (2-dimensional 2-layer network) -------------------------- #define I_CX 1 // 0 - No layer, 1 - Add Layer of Pyramidal neurons #define I_IN 1 // 0 - No layer, 1 - Add Layer of Interneurons #define Mcx 256 // Number of CX cells (Pyramidal neurons) in X direction #define Min 128 // Number of IN cells (Interneurons) in X direction #define Mcx1 256 // Number of CX cells in Y direction #define Min1 128 // Number of IN cells in Y direction #define Max ((M > Mcx) ? M : Mcx) #define Max1 ((M1 > Mcx1) ? M1 : Mcx1) #define SELFINGcx 0 // 0 - CX without self excitation; 1 - with ... //-------------- Define the connections between cells ----------------- #define MS_CX_CX 8 //radius of CX->CX connections in X direction #define MS_CX_CX1 8 //radius of CX->CX connections in Y direction #define MS_CX_CX_MAX ((MS_CX_CX > MS_CX_CX1) ? MS_CX_CX : MS_CX_CX1) #define N_CX_CX (2*MS_CX_CX+1)*(2*MS_CX_CX1+1) //total number of CX connections accepted by one CX #define MS_CX_IN 8 //radius of CX->IN connections in X direction #define MS_CX_IN1 8 //radius of CX->IN connections in Y direction #define MS_CX_IN_MAX ((MS_CX_IN > MS_CX_IN1) ? MS_CX_IN : MS_CX_IN1) #define N_CX_IN (2*MS_CX_IN+1)*(2*MS_CX_IN1+1) //total number of CX connections accepted by one IN #define MS_IN_CX 2 //radius of IN->CX connections in X direction #define MS_IN_CX1 2 //radius of IN->CX connections in Y direction #define MS_IN_CX_MAX ((MS_IN_CX > MS_IN_CX1) ? MS_IN_CX : MS_IN_CX1) #define N_IN_CX (2*MS_IN_CX+1)*(2*MS_IN_CX1+1) //-------------------RS CELL (main class to describe regular spiking neuron)----------------------- class RS{ double xp, xpp, mu, sigma_n, beta_n; double sigma_dc, beta_dc; public: double x, y, alpha, sigma, sigma_e, beta_e, Idc, S_CX_DEND; int spike; RS(){ mu=0.0005; spike=0; alpha=3.65; sigma=6.00E-2; sigma_e=1.0; sigma_dc=1.0; beta_e=0.133; beta_dc=0.133; Idc = 0; S_CX_DEND = 165.0e-6; xpp=-1+sigma; xp=-1+sigma; x=-1+sigma; y= x - alpha/(1-x); } void init(){ mu=0.0005; spike=0; alpha=3.65; sigma=6.00E-2; sigma_e=1.0; sigma_dc=1.0; beta_e=0.133; beta_dc=0.133; Idc = 0; S_CX_DEND = 165.0e-6; xpp=-1+sigma; xp=-1+sigma; x=-1+sigma; y= x - alpha/(1-x); } void calc(double); }; void RS::calc(double I){ beta_n = beta_e * I + beta_dc * Idc; sigma_n = sigma_e * I + sigma_dc * Idc; if (beta_n < -1.0) beta_n = -1.0; if (beta_n > 1.0) beta_n = 1.0; if(xp <= 0.0) { x = alpha / (1.0 - xp) + y +beta_n; spike = 0; } else{ if(xp <= alpha + y +beta_n && xpp <= 0.0) { x = alpha + y + beta_n; spike = 1; } else { x = -1; spike = 0; } } y = y - mu* (xp +1.0) + mu * sigma + mu * sigma_n; xpp = xp; xp = x; // y=1; } //-------------------FS1 CELL (main class to describe fast spiking interneuron)----------------- class FS1{ double xp, xpp, mu, sigma_n, beta_n, ii, gg; double sigma_dc, beta_dc, beta_ii, sigma_ii; public: double x, y, alpha, sigma, sigma_e, beta_e, Idc, S_CX_DEND; int spike; FS1(){ ii=0; mu=0.002; spike=0; alpha=3.8; //3.87; sigma=-1.75E-2; sigma_e=1.0; sigma_dc=1.0; beta_e=0.1; beta_dc=0.0; Idc = 0; S_CX_DEND = 165.0e-6; xpp=-1+sigma; xp=-1+sigma; x=-1+sigma; y= x - alpha/(1-x); gg = 0.5; beta_ii = 0.5; sigma_ii = 0.5; } void init(){ ii=0; mu=0.002; spike=0; alpha=3.8; //3.87; sigma=-1.75E-2; sigma_e=1.0; sigma_dc=1.0; beta_e=0.1; beta_dc=0.0; Idc = 0; S_CX_DEND = 165.0e-6; xpp=-1+sigma; xp=-1+sigma; x=-1+sigma; y= x - alpha/(1-x); gg = 0.5; beta_ii = 0.5; sigma_ii = 0.5; } void calc(double); }; void FS1::calc(double I){ if(spike > 0.1){ ii = 0.60 * ii - gg; } else{ ii = 0.60 * ii; } beta_n = beta_e * I + beta_dc * Idc; sigma_n = sigma_e * I + sigma_dc * Idc; if (beta_n < -0.0001) beta_n = -0.0001; if (beta_n > 1.0) beta_n = 1.0; beta_n = beta_n + beta_ii *ii; sigma_n = sigma_n + sigma_ii *ii; if(xp <= 0.0) { x = alpha / (1.0 - xp) + y +beta_n; spike = 0; } else{ if(xp <= alpha + y +beta_n && xpp <= 0.0) { x = alpha + y + beta_n; spike = 1; } else { x = -1; spike = 0; } } // y = y - mu* (xp +1.0) + mu * sigma + mu * sigma_n; y = -2.9; xpp = xp; xp = x; } //------------Main class to describe first order kinetic model for AMPA map-based------------ //--------------------------------------synapse with depression-------------------------------- class AMPAmapD { static double E_AMPA; static double gamma; double d; // depression variable public: double I, d_dep, d_rec; AMPAmapD() { I=0; d=1; d_dep=0.5; d_rec=0.010; } void calc(double, double, int); }; double AMPAmapD::E_AMPA = 0, AMPAmapD::gamma = 0.6; void AMPAmapD::calc(double g_AMPA, double x_post, int spike) { if(spike > 0.1){ I = gamma * I - d*g_AMPA * (x_post - E_AMPA); d = (1.0 - d_dep)*d; } else{ I = gamma * I; d = 1.0 - (1.0 - d_rec) * (1.0 - d); } } //---------Main class to describe first order kinetic model for AMPA map-based synapse------- class AMPAmap { static double E_AMPA; static double gamma; public: double I; AMPAmap() { I=0; } void calc(double, double, int); }; double AMPAmap::E_AMPA = 0, AMPAmap::gamma = 0.6; void AMPAmap::calc(double g_AMPA, double x_post, int spike) { if(spike > 0.1){ I = gamma * I - g_AMPA * (x_post - E_AMPA); } else{ I = gamma * I; } } //--------Main class to describe first order kinetic model for GABA-A map-based synapse---------- class GABAAmap { static double E_GABAA; static double gamma; public: double I; GABAAmap() { I=0; } void calc(double, double, int); }; double GABAAmap::E_GABAA = -1.1, GABAAmap::gamma = 0.8; void GABAAmap::calc(double g_GABA, double x_post, int spike) { if(spike > 0.1){ I = gamma * I - g_GABA * (x_post - E_GABAA); } else{ I = gamma * I; } } //-----Main class to describe first order kinetic model for AMPA synapse used for external--------- //--------------------------- stimulation---------------------------------------------------------- class Extern_ampa { static double Cdur, Cmax, Deadtime, Prethresh; double R, C, R0, R1; double lastrelease; double q, Rinf, Rtau; double TR, w, wom, RRR; double exptable(double z) { if((z > -10) && (z < 10)) return( exp(z) ); else return( 0 ); } public: double g, Alpha, Beta; Extern_ampa() { Alpha = 0.94; Beta = 0.18; R = 0, C = 0, R0 = 0, R1 = 0; lastrelease = -100; Rinf = Cmax*Alpha / (Cmax*Alpha + Beta); Rtau = 1 / ((Alpha * Cmax) + Beta); TR = 1000, w=0.01, wom=0; } void iii(unsigned int seek) {srand(seek);} void calc(double g_Extern_ampa, double x); }; double Extern_ampa::Cdur = 0.3, Extern_ampa::Cmax = 0.5, Extern_ampa::Deadtime = 1; double Extern_ampa::Prethresh = 0; void Extern_ampa::calc(double g_Extern_ampa, double x) { q = ((x - lastrelease) - Cdur); if (q > Deadtime) { if ((x - lastrelease) > TR) { C = Cmax; R0 = R; lastrelease = x; // RRR = 1.0*rand()/(RAND_MAX + 1.0); // if(RRR < 0.000001) RRR = 0.000001; // TR = -(log(RRR))/(w+(w/2)*sin(wom*x)); } } else if (q < 0) { } else if (C == Cmax) { R1 = R; C = 0.; } if (C > 0) { R = Rinf + (R0 - Rinf) * exptable (-(x - lastrelease) / Rtau); } else { R = R1 * exptable (-Beta * (x - (lastrelease + Cdur))); } g = g_Extern_ampa * R; } //+++++++++++++++++++ MAIN PROGRAM +++++++++++++++++++++++++++++++++++++++++++ //----------external functions---------------------------------------------- void fun(double); //----------external variables --------------------------------------------- double *g_ampa_cx_cx[Mcx][Mcx1], *g_ampa_cx_in[Min][Min1], *g_gaba_a_in_cx[Mcx][Mcx1]; double g_ext_cx, g_ext_in; FILE *f1, *f11, *f13, *f98, *f99; int C_CXCX[2*MS_CX_CX+1][2*MS_CX_CX1+1]; int C_CXIN[2*MS_CX_IN+1][2*MS_CX_IN1+1]; int C_INCX[2*MS_IN_CX+1][2*MS_IN_CX1+1]; int k_CXCX, k_CXIN, k_INCX; int k_CXCXmax=0, k_CXINmax=0, k_INCXmax=0; //----------external classes (beginning of initialization)------------------ //---Synapses--------------------- AMPAmapD *a_cx_cx[Mcx][Mcx1]; AMPAmap *a_cx_in[Min][Min1]; GABAAmap *ga_in_cx[Mcx][Mcx1]; //---Neurons----------------------- RS cx_cell[Mcx][Mcx1]; FS1 in_cell[Min][Min1]; //---External input---------------- Extern_ampa a_ext1; //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ main(int argc,char **argv) { //---------general parameters---------------------------------------------- double t = 0, tmax, t3D, ttime; double h = 0.02, R, av=0; int i, j, i1, j1, k, ii = 0, jj=0, i_gip = 0; double scale; double g_AMPA_CX_CX, g_AMPA_CX_IN, g_GABA_A_IN_CX; double g_Extern_ampa; double TTR=0, aver=0; int TTRS=0; char namein[10], namecx[10]; //-------parameter initialization (from file)---------------------------------- f1=fopen("input2D.txt","r"); fscanf(f1, "%lf %lf %lf %lf %lf %lf %lf ", &tmax, &t3D, &ttime, &g_AMPA_CX_CX, &g_AMPA_CX_IN, &g_GABA_A_IN_CX, &g_Extern_ampa); printf("\n param: %lf %lf %lf AM_CX_CX=%lf AM_CX_IN=%lf GA_IN_CX=%lf %lf ", tmax, t3D, ttime, g_AMPA_CX_CX, g_AMPA_CX_IN, g_GABA_A_IN_CX, g_Extern_ampa); fclose(f1); //----------classes initialization (continue)---------------------------- printf("\n CLASS INIT CONT"); for(i=0; i < Mcx; ++i) for(j=0; j < Mcx1; ++j){ a_cx_cx[i][j] = new AMPAmapD[N_CX_CX]; ga_in_cx[i][j] = new GABAAmap[N_IN_CX]; g_ampa_cx_cx[i][j] = new double[N_CX_CX]; g_gaba_a_in_cx[i][j] = new double[N_IN_CX]; } for(i=0; i < Min; ++i) for(j=0; j < Min1; ++j){ a_cx_in[i][j] = new AMPAmap[N_CX_IN]; g_ampa_cx_in[i][j] = new double[N_CX_IN]; } //--------set synaptic conductances to zero (will be updated later)---------------- if(I_CX == 1){ for(i=0; i < Mcx; ++i) for(j=0; j < Mcx1; ++j) for(k=0; k < N_CX_CX; ++k) g_ampa_cx_cx[i][j][k] = 0; for(i=0; i < Mcx; ++i) for(j=0; j < Mcx1; ++j) for(k=0; k < N_IN_CX; ++k) g_gaba_a_in_cx[i][j][k] = 0; g_ext_cx = 0; } if(I_IN == 1){ for(i=0; i < Min; ++i) for(j=0; j < Min1; ++j) for(k=0; k < N_CX_IN; ++k) g_ampa_cx_in[i][j][k] = 0; g_ext_in = 0; } //----------here we are changing some variables to introduce random variability--------- srand(6); for(i=0; i < Mcx; ++i) for(j=0; j < Mcx1; ++j){ R = 2.0 * rand()/(RAND_MAX + 1.0) - 1.0; cx_cell[i][j].sigma = cx_cell[i][j].sigma + R * 0.02; } cout << "\n -------------------------------------------------"; printf("\n VARIABILITY"); printf("\n VARIABILITY"); //--------------open ALL files------------------------------------- f11 = fopen("time_cx", "w"); f13 = fopen("time_in", "w"); sprintf(namecx,"graf_cx%d", TTRS); f98 = fopen(namecx, "w"); sprintf(namein,"graf_in%d", TTRS); f99 = fopen(namein, "w"); //----------Connection matrix------------------------------- printf("\n Begin Connect Matrix"); //---set connection matrix to zero-------------------------- for(i=0; i<(2*MS_CX_CX+1); i++) for(j=0; j<(2*MS_CX_CX1+1); j++){ k_CXCX=0; C_CXCX[i][j]=0; } for(i=0; i<(2*MS_CX_IN+1); i++) for(j=0; j<(2*MS_CX_IN1+1); j++){ k_CXIN=0; C_CXIN[i][j]=0; } for(i=0; i<(2*MS_IN_CX+1); i++) for(j=0; j<(2*MS_IN_CX1+1); j++){ k_INCX=0; C_INCX[i][j]=0; } //---start updating connection matrix------------------------- for(i1=-MS_CX_CX; i1<=MS_CX_CX; i1++) for(j1=-MS_CX_CX1; j1<=MS_CX_CX1; j1++){ scale = sqrt((double) (i1*i1 + j1*j1)); if(scale <= MS_CX_CX_MAX){ C_CXCX[i1+MS_CX_CX][j1+MS_CX_CX1]=1; k_CXCX=k_CXCX+1;} if( (scale == 0) && (SELFINGcx == 0) ){ C_CXCX[i1+MS_CX_CX][j1+MS_CX_CX1]=0; k_CXCX=k_CXCX-1;} } for(i1=-MS_IN_CX; i1<=MS_IN_CX; i1++) for(j1=-MS_IN_CX1; j1<=MS_IN_CX1; j1++){ scale = sqrt((double) (i1*i1 + j1*j1)); if(scale <= MS_IN_CX_MAX){ C_INCX[i1+MS_IN_CX][j1+MS_IN_CX1]=1; k_INCX=k_INCX+1;} } for(i1=-MS_CX_IN; i1<=MS_CX_IN; i1++) for(j1=-MS_CX_IN1; j1<=MS_CX_IN1; j1++){ scale = sqrt((double) (i1*i1 + j1*j1)); if(scale <= MS_CX_IN_MAX){ C_CXIN[i1+MS_CX_IN][j1+MS_CX_IN1]=1; k_CXIN=k_CXIN+1;} } k_CXCXmax=k_CXCX; k_CXINmax=k_CXIN; k_INCXmax=k_INCX; printf("\n End Connect Matrix"); //----------------CALCULATION---------------------------------------- printf("\n CALCULATION IN PROGRESS!!!: t= %lf: tmax= %lf", t,tmax); t=0; printf("\n start simulation"); while( t < tmax){ t++; fun(t); // External (DC type) stimulation to 10 pyramidal neurons (only applied once) for(i=0; i<10; i++){ if(t>300) cx_cell[i][Mcx1-1].Idc=0.2; if(t>320) cx_cell[i][Mcx1-1].Idc=0; } //------- Calculate average response-------------- aver=0; for(i=Mcx/2-4; i<=Mcx/2+5; i++) for(j=Mcx1/2-4; j<=Mcx1/2+5; j++){ aver=aver+cx_cell[i][j].x; } //--------scaling of all conductances----------------------------------- if( (t>200) && (t<=201.5) ){ printf("\n CONDUCTANCE INITIALIZATION"); if(I_CX == 1){ g_ampa_cx_cx[0][0][0] = g_AMPA_CX_CX / cx_cell[0][0].S_CX_DEND; g_gaba_a_in_cx[0][0][0] = g_GABA_A_IN_CX / cx_cell[0][0].S_CX_DEND; } if(I_IN == 1){ g_ampa_cx_in[0][0][0] = g_AMPA_CX_IN / in_cell[0][0].S_CX_DEND; } printf("\n COUPLING %lf %lf %lf",g_ampa_cx_cx[0][0][0], g_gaba_a_in_cx[0][0][0], g_ampa_cx_in[0][0][0]); } //----------------------------------------------------------------------- av=0; if((t/(5000))*(5000) == t){ printf("\n %lf", t/tmax); } //----Print output to files---------------- if((t > ttime) && (((int)t/(2))*(2) == (int)t)){ if(I_CX == 1) fprintf(f11,"%lf ",t); if(I_IN == 1) fprintf(f13,"%lf ",t); for(i = 0; i < Mcx; ++i){ if(I_CX == 1) fprintf(f11,"%lf ", cx_cell[i][0].x); } for(i = 0; i < Min; ++i){ if(I_IN == 1) fprintf(f13,"%lf ", in_cell[i][0].x); } fprintf(f11,"\n"); fprintf(f13,"\n"); } if((t > t3D) && (((int)t/(10))*(10) == (int)t)){ for(i = 0; i < Mcx; ++i){ for(j = 0; j < Mcx1; ++j){ if(I_CX == 1) fprintf(f98,"%lf ", cx_cell[i][j].x); } fprintf(f98,"\n"); } for(i = 0; i < Min; ++i){ for(j = 0; j < Min1; ++j){ if(I_IN == 1) fprintf(f99,"%lf ", in_cell[i][j].x); } fprintf(f99,"\n"); } } } //--------------------END CALCULATION------------------------------- //-----------------close ALL files----------------------------------- fclose(f11); fclose(f13); fclose(f98); fclose(f99); printf("\n"); } //+++++++++++ Function to calculate the right sides for ALL ODE +++++++++++++++ void fun(double t){ int i, j, k, kmax, i1, j1, ii, jj; double IsynCXCX[Mcx][Mcx1], IsynCXIN[Min][Min1], IsynINCX[Mcx][Mcx1]; for(i = 0; i < Mcx; ++i) for(j = 0; j < Mcx1; ++j){ //-----calculating reciprocal AMPA-mediated synaptic currents between CX cells-------------- if(I_CX == 1){ for(i1 = -MS_CX_CX, k = 0; i1 <= MS_CX_CX; ++i1) for(j1 = -MS_CX_CX1; j1 <= MS_CX_CX1; ++j1){ ii=i+i1, jj=j+j1; if((ii>=0) && (ii=0) && (jj 0){ a_cx_cx[i][j][k].calc(g_ampa_cx_cx[0][0][0], cx_cell[i][j].x, cx_cell[ii][jj].spike); ++k; } }} kmax = k; IsynCXCX[i][j] = 0; for(k = 0; k < kmax; ++k){ if(kmax == 0) { } //NO synapses else { IsynCXCX[i][j] = IsynCXCX[i][j] + a_cx_cx[i][j][k].I / k_CXCXmax;} } } //------calculating GABA-A type current from IN to CX cells-------------------------------- if((I_CX == 1) && (I_IN == 1)){ for(i1 = -MS_IN_CX, k = 0; i1 <= MS_IN_CX; ++i1) for(j1 = -MS_IN_CX1; j1 <= MS_IN_CX1; ++j1){ ii=i*Min/Mcx+i1, jj=j*Min1/Mcx1+j1; if((ii>=0) && (ii=0) && (jj 0){ ga_in_cx[i][j][k].calc(g_gaba_a_in_cx[0][0][0], cx_cell[i][j].x, in_cell[ii][jj].spike); ++k; } } } kmax = k; IsynINCX[i][j] = 0; for(k = 0; k < kmax; ++k){ if(kmax == 0) { } //NO synapses else { IsynINCX[i][j] = IsynINCX[i][j] + ga_in_cx[i][j][k].I / k_INCXmax;} } } } for(i = 0; i < Min; ++i) for(j = 0; j < Min1; ++j){ //------calculating AMPA type current from CX to IN cells-------------------------------- if((I_CX == 1) && (I_IN == 1)){ for(i1 = -MS_CX_IN, k = 0; i1 <= MS_CX_IN; ++i1) for(j1 = -MS_CX_IN1; j1 <= MS_CX_IN1; ++j1){ ii=i*Mcx/Min+i1, jj=j*Mcx1/Min1+j1; if((ii>=0) && (ii=0) && (jj 0){ a_cx_in[i][j][k].calc(g_ampa_cx_in[0][0][0], in_cell[i][j].x, cx_cell[ii][jj].spike); ++k; } } } kmax = k; IsynCXIN[i][j] = 0; for(k = 0; k < kmax; ++k){ if(kmax == 0) { } //NO synapses else { IsynCXIN[i][j] = IsynCXIN[i][j] + a_cx_in[i][j][k].I / k_CXINmax; } } } } //------calculating the state of all neurons at the next iteration based on synaptic input----- for(i=0; i < Mcx; ++i) for(j=0; j < Mcx1; ++j){ if(I_CX == 1) cx_cell[i][j].calc(IsynCXCX[i][j]+IsynINCX[i][j]); } for(i=0; i < Min; ++i) for(j=0; j < Min1; ++j){ if(I_IN == 1) in_cell[i][j].calc(IsynCXIN[i][j]); } }