Large cortex model with map-based neurons (Rulkov et al 2004)

 Download zip file 
Help downloading and running models
Accession:45525
We develop a new computationally efficient approach for the analysis of complex large-scale neurobiological networks. Its key element is the use of a new phenomenological model of a neuron capable of replicating important spike pattern characteristics and designed in the form of a system of difference equations (a map). ... Interconnected with synaptic currents these model neurons demonstrated responses very similar to those found with Hodgkin-Huxley models and in experiments. We illustrate the efficacy of this approach in simulations of one- and two-dimensional cortical network models consisting of regular spiking neurons and fast spiking interneurons to model sleep and activated states of the thalamocortical system. See paper for more.
Reference:
1 . Rulkov NF, Timofeev I, Bazhenov M (2004) Oscillations in large-scale cortical networks: map-based model. J Comput Neurosci 17:203-23 [PubMed]
Citations  Citation Browser
Model Information (Click on a link to find other models with that property)
Model Type: Realistic Network;
Brain Region(s)/Organism: Neocortex;
Cell Type(s): Neocortex L5/6 pyramidal GLU cell; Neocortex fast spiking (FS) interneuron;
Channel(s):
Gap Junctions:
Receptor(s): GabaA; AMPA;
Gene(s):
Transmitter(s):
Simulation Environment: C or C++ program;
Model Concept(s): Activity Patterns; Oscillations; Spatio-temporal Activity Patterns; Simplified Models; Sleep;
Implementer(s): Bazhenov, Maxim [Bazhenov at Salk.edu]; Rulkov, Nikolai [nrulkov at ucsd.edu];
Search NeuronDB for information about:  Neocortex L5/6 pyramidal GLU cell; GabaA; AMPA;
/
largectx2004
index.html
colors.m
import2D.m
input2D.txt
movieMAP.m
network2D.cpp
                            
//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 <stdio.h> 
#include <stdlib.h> 
#include <math.h> 
#include <iostream.h> 

//---------------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 < -0.0001) beta_n = -0.0001;
     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.1) && (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 / k_CXCXmax;
  g_gaba_a_in_cx[0][0][0] = g_GABA_A_IN_CX / cx_cell[0][0].S_CX_DEND / k_INCXmax;
}
if(I_IN == 1){
  g_ampa_cx_in[0][0][0] = g_AMPA_CX_IN / in_cell[0][0].S_CX_DEND / k_CXINmax;
}
 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<Mcx) && (jj>=0) && (jj<Mcx1)){ 
       if(C_CXCX[i1+MS_CX_CX][j1+MS_CX_CX1] > 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;
   } 
  } 
 } 

//------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<Min) && (jj>=0) && (jj<Min1)){ 
              if(C_INCX[i1+MS_IN_CX][j1+MS_IN_CX1] > 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;
  } 
  } 
} 
   
 }
 
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<Mcx) && (jj>=0) && (jj<Mcx1)){ 
            if(C_CXIN[i1+MS_CX_IN][j1+MS_CX_IN1] > 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;

   } 
  } 
 }   
 }

//------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]); 
  } 

}