A model of antennal lobe of bee (Chen JY et al. 2015)

 Download zip file 
Help downloading and running models
Accession:226471
" ... Here we use calcium imaging to reveal how responses across antennal lobe projection neurons change after association of an input odor with appetitive reinforcement. After appetitive conditioning to 1-hexanol, the representation of an odor mixture containing 1-hexanol becomes more similar to this odor and less similar to the background odor acetophenone. We then apply computational modeling to investigate how changes in synaptic connectivity can account for the observed plasticity. Our study suggests that experience-dependent modulation of inhibitory interactions in the antennal lobe aids perception of salient odor components mixed with behaviorally irrelevant background odors."
Reference:
1 . Chen JY, Marachlian E, Assisi C, Huerta R, Smith BH, Locatelli F, Bazhenov M (2015) Learning modifies odor mixture processing to improve detection of relevant components. J Neurosci 35:179-97 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Realistic Network;
Brain Region(s)/Organism:
Cell Type(s): Hodgkin-Huxley neuron;
Channel(s):
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: C or C++ program;
Model Concept(s): Learning; Olfaction;
Implementer(s): Chen, Jen-Yung [chen.jenyung at gmail.com]; Assisi, Collins ;
#include <stdio.h>
#include <stdlib.h> 
#include <math.h>
#include <iostream>
using namespace std;

//---------------CONSTANTs initialization----------------------------
//--------------- Network Geometry ------------------------------------
#define I_TC    1     //  0 - No layer, 1 - Add Layer
#define I_RE    1     //  0 - No layer, 1 - Add Layer
#define I_CX    0     //  0 - No layer, 1 - Add Layer
#define I_IN    0     //  0 - No layer, 1 - Add Layer
#define I_GB    1     //  0 - No GABAb from RE to TC, 1 - Yes
#define I_GB_RERE    0     //  0 - No GABAb from RE to RE, 1 - Yes

#define Mre       280      //  Number of LN in olfaction
#define Mtc       100      //  Number of PN in olfaction
#define Mc      1      //  Number of CX-IN cells in X direction
#define Mre1      1      //  Number of RE cells in Y direction
#define Mtc1      1      //  Number of TC cells in Y direction
#define Mc1     1      //  Number of CX-IN cells in Y direction
#define Max     ((Mre > Mtc) ? Mre : Mtc)
#define Max1    ((Mre1 > Mtc1) ? Mre1 : Mtc1)

//-------------Boundary conditions ---------------------------------------
#define BOUND    0   //  0 - flow; 1 - periodic; 9 - no boundary elements 
#define SELFING    0   //  0 - RE without self inhibition; 1 - with ... 
#define SELFINGcx  0   //  0 - CX without self excitation; 1 - with ... 

//-------------- Define the connections between cells -----------------
#define MS_RE_RE 0
#define MS_RE_RE1 0
#define MS_RE_RE_MAX  ((MS_RE_RE > MS_RE_RE1) ? MS_RE_RE : MS_RE_RE1)
#define N_RE_RE  Mre*Mre1 //(2*MS_RE_RE+1)*(2*MS_RE_RE1+1) 

#define MS_RE_TC 3
#define MS_RE_TC1 0
#define MS_RE_TC_MAX ((MS_RE_TC > MS_RE_TC1) ? MS_RE_TC : MS_RE_TC1)
#define N_RE_TC  Mre*Mre1 //(2*MS_RE_TC+1)*(2*MS_RE_TC1+1) //number of
                          //connections accepted by some TC from all RE

#define MS_TC_RE 3
#define MS_TC_RE1 0
#define MS_TC_RE_MAX ((MS_TC_RE > MS_TC_RE1) ? MS_TC_RE : MS_TC_RE1)
#define N_TC_RE  Mtc*Mtc1 //(2*MS_TC_RE+1)*(2*MS_TC_RE1+1)

#define MS_TC_TC 3
#define MS_TC_TC1 0
#define MS_TC_TC_MAX ((MS_TC_TC > MS_TC_TC1) ? MS_TC_TC : MS_TC_TC1)
#define N_TC_TC  Mtc*Mtc1 //(2*MS_TC_TC+1)*(2*MS_TC_TC1+1)

#define MS_CX_CX 0
#define MS_CX_CX1 0
#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)

#define MS_CX_IN 0 
#define MS_CX_IN1 0
#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)

#define MS_IN_CX 0
#define MS_IN_CX1 0
#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)

#define MS_TC_CX 0
#define MS_TC_CX1 0
#define MS_TC_CX_MAX  ((MS_TC_CX > MS_TC_CX1) ? MS_TC_CX : MS_TC_CX1)
#define N_TC_CX  (2*MS_TC_CX+1)*(2*MS_TC_CX1+1)

#define MS_TC_IN 0
#define MS_TC_IN1 0
#define MS_TC_IN_MAX  ((MS_TC_IN > MS_TC_IN1) ? MS_TC_IN : MS_TC_IN1)
#define N_TC_IN  (2*MS_TC_IN+1)*(2*MS_TC_IN1+1)

#define MS_CX_TC 0 
#define MS_CX_TC1 0
#define MS_CX_TC_MAX  ((MS_CX_TC > MS_CX_TC1) ? MS_CX_TC : MS_CX_TC1)
#define N_CX_TC  (2*MS_CX_TC+1)*(2*MS_CX_TC1+1)

#define MS_CX_RE 0 
#define MS_CX_RE1 0
#define MS_CX_RE_MAX  ((MS_CX_RE > MS_CX_RE1) ? MS_CX_RE : MS_CX_RE1)
#define N_CX_RE  (2*MS_CX_RE+1)*(2*MS_CX_RE1+1)

//------------Number of ODE for earch cell -------------------------------
#define N_RE 9 //4 //7
#define N_TC 6 //12 
#define N_GB 2
#define N_GA 1

#define N_DEND   8
#define N_SOMA   3
#define N_CX     (N_DEND + N_SOMA)
#define N_IN     N_CX   //4

#define N_EQ1  (N_RE*I_RE + N_RE_RE*N_GB*I_RE*I_GB_RERE)*Mre*Mre1
#define N_EQ2  (N_TC*I_TC + N_RE_TC*N_GB*I_RE*I_TC*I_GB)*Mtc*Mtc1
#define N_EQ3  (N_CX*I_CX + N_IN*I_IN + N_IN_CX*N_GB*I_CX*I_IN)*Mc*Mc1
#define N_EQ4  (N_RE_TC*N_GA*I_RE*I_TC)*Mtc*Mtc1
#define N_EQ5  (N_RE_RE*N_GA*I_RE)*Mre*Mre1
#define N_EQ   (N_EQ1+N_EQ2+N_EQ3+N_EQ4+N_EQ5)    //  Complete number of ODE

//++++++++++++++ CURRENTs DESCRIPTION ++++++++++++++++++++++++++++++++++++++
//---------------Low-threshold Ca2+ current (RE cell)---------------------
class IT_RE {
  static double Shift, Ca_0, Cels;
  double m_inf, tau_m, h_inf, tau_h, ratio, eca, Phi_m, Phi_h, eca0;
public:
  double iT, m0, h0, Qm, Qh;
  double G_Ca;
  IT_RE(double v) {
    G_Ca = 1.75;
    Qm = 2.5; //2.5; //3; 
    Qh = 3; //2.5; //5;
    Phi_m = pow(Qm,((Cels-24)/10));
    Phi_h = pow(Qh,((Cels-24)/10));
    m0 = 1/(1 + exp(-(v + Shift + 50)/7.4));
    h0 = 1/(1 + exp((v + Shift + 78)/5));
    eca0 = 1000*8.31441*(273.15 + Cels)/(2*96489);
    } 
  void calc(double m, double h, double &fm, double &fh, 
            double v, double cai, double x);
};

double IT_RE::Shift = 2, IT_RE::Ca_0 = 2, IT_RE::Cels = 36;

void IT_RE::calc(double m, double h, double &fm, double &fh,
                 double v, double cai, double x) {
  ratio = Ca_0/cai;
    if(ratio <= 0.)printf("\n LOG ERROR: RE: cai=%lf ratio=%lf",cai,ratio);
  eca = eca0 * log(ratio);
  iT = G_Ca*m*m*h*(v - eca);                              
  m_inf = 1/(1 + exp(-(v + 52)/7.4));
  tau_m = (3 + 1/(exp((v + 27)/10) + exp(-(v + 102)/15)))/Phi_m;
  h_inf = 1/(1 + exp((v + 80)/5));
  tau_h = (85 + 1/(exp((v + 48)/4) + exp(-(v + 407)/50)))/Phi_h;
  fm = -(m - m_inf)/tau_m;                                  
  fh = -(h - h_inf)/tau_h;                                  
}

//-------------- Ca-dependent potassium current (RE cell) --------------------
class IKCa_RE {
  static double E_KCa, Alpha, Beta, Cels;     
  double m_inf, tau_m, Tad, car;                                 
public:
  double iKCa, m0;
  double G_KCa;
  IKCa_RE(double cai) {
    G_KCa = 0.; //7; 
    Tad = pow(3,((Cels-22)/10));
    car = Alpha*cai*cai/Beta;
    m0 = car/(car + 1);  } 
  void calc(double m, double &fm, double v, double cai, double x);
};

double IKCa_RE::E_KCa = -95;
double IKCa_RE::Alpha = 48, IKCa_RE::Beta = 0.03, IKCa_RE::Cels = 36;   

void IKCa_RE::calc(double m, double &fm, double v, double cai, double x){
  iKCa = G_KCa*m*m*(v - E_KCa);                         
  car = Alpha*cai*cai/Beta;
  m_inf = car/(car + 1);
  tau_m = (1/(Beta*(car + 1)))/Tad;
  if(tau_m < 0.1) tau_m = 0.1;
  fm = -(1/tau_m)*(m - m_inf);                   
}

//-------------- Non-specific current (RE cell) --------------------
class ICAN_RE {
  static double E_CAN, Alpha, Beta, Cels;   
  double m_inf, tau_m, Tad, car;                                 
public:
  double iCAN, m0;
  double G_CAN;
  ICAN_RE(double cai) {
    G_CAN = 0; //0.06; 
    Tad = pow(3,((Cels-22)/10));
    car = Alpha*cai*cai/Beta;
    m0 = car/(car + 1);    } 
  void calc(double m, double &fm, double v, double cai, double x);
};

double ICAN_RE::E_CAN = -20;
double ICAN_RE::Alpha = 20, ICAN_RE::Beta = 0.002, ICAN_RE::Cels = 36;          

void ICAN_RE::calc(double m, double &fm, double v, double cai, double x){
  iCAN = G_CAN*m*m*(v - E_CAN);  
  car = Alpha*cai*cai/Beta;
  m_inf = car/(car + 1);
  tau_m = (1/(Beta*(car + 1)))/Tad;
  if(tau_m  < 0.1) tau_m = 0.1;
  fm = -(1/tau_m)*(m - m_inf); 
}

//--------------fast Na and K current (RE and TC cells)------------------
class INaK {
  static double Cels;
  double Alpha1, Beta1, Alpha2, Beta2, Alpha3, Beta3, v2, v2K, Phi;
  double tau_m, m_inf, tau_h, h_inf, tau_n, n_inf;
public:
  static double E_Na, E_K;
  double iK, iNa, m0, h0, n0;
  double G_Na, G_K, Vtr, VtrK;
  double S1, S2;
  INaK(double v) {
    G_K = 10;///////////////////////
    G_Na = 100;/////////////////////
    Vtr = -50;
    VtrK = -50;
    S1 = 0.32;
    S2 = 0.02;
    v2 = v - Vtr;
    v2K = v - VtrK;
    Phi = pow(3,((Cels-36)/10));
    Alpha1 = 0.32*(13 - v2)/(exp((13 - v2)/4) - 1);
    Beta1 = 0.28*(v2 - 40)/(exp((v2 - 40)/5) - 1);
    m0 = Alpha1/(Alpha1 + Beta1);

    Alpha2 = 0.128*exp((17 - v2)/18);
    Beta2 = 4/(exp((40 - v2)/5) + 1);
    h0 = Alpha2/(Alpha2 + Beta2);

    Alpha3 = 0.02*(15 - v2)/(exp((15 - v2)/5) - 1);
    Beta3 = 0.5*exp((10 - v2)/40);
    n0 = Alpha3/(Alpha3 + Beta3);     } 
  void calc(double m, double h, double n, double &fm, double &fh, double &fn, 
            double v, double x);
};

double INaK::E_K = -95, INaK::E_Na = 50, INaK::Cels = 22; 

void INaK::calc(double m, double h, double n, double &fm, double &fh, double &fn,
                   double v, double x){
  v2 = v - Vtr;
  v2K = v - VtrK;
  iNa = G_Na*m*m*m*h*(v - E_Na);
  Alpha1 = S1*(13 - v2)/(exp((13 - v2)/4) - 1);
  Beta1 = 0.28*(v2 - 40)/(exp((v2 - 40)/5) - 1);
  tau_m = 1/(Alpha1 + Beta1) / Phi;
  m_inf = Alpha1/(Alpha1 + Beta1);

  Alpha2 = 0.128*exp((17 - v2)/18);
  Beta2 = 4/(exp((40 - v2)/5) + 1);
  tau_h = 1/(Alpha2 + Beta2) / Phi;
  h_inf = Alpha2/(Alpha2 + Beta2);

  fm = -(m - m_inf)/tau_m;                 
  fh = -(h - h_inf)/tau_h;                 

  iK = G_K* n*n*n*n*(v - E_K);    
  Alpha3 = S2*(15 - v2K)/(exp((15 - v2K)/5) - 1);
  Beta3 = 0.5*exp((10 - v2K)/40);
  tau_n = 1/(Alpha3 + Beta3) / Phi;
  n_inf = Alpha3/(Alpha3 + Beta3);
  
  fn  = -(n - n_inf)/tau_n;                 
}

//------------------Ca-dynamics------------------------------------
class ICa {
  static double Ca_inf, K_T, K_d;
  double drive, drive0;                                 
public:
  double Taur, D;
  ICa() {Taur = 5; D = 0.1;
         drive0 = 10.0/(2.*96489.); }
  void calc(double cai, double &fcai, double iT, double x);
};

double ICa::Ca_inf = 2.4e-4;
double ICa::K_T = 0.0001, ICa::K_d = 0.0001;

void ICa::calc(double cai, double &fcai, double iT, double x) {
  drive = -drive0 * iT / D;
  if(drive < 0) drive = 0;
  fcai = drive + (Ca_inf - cai)/Taur; // - K_T*cai/(cai + K_d);
}

//------------------Low-theshold Ca2+ current (TC cell)-----------------
class IT_TC {
  static double Ca_0, Cels, Qm, Qh, Shift;
  double m_inf, tau_m, h_inf, tau_h, Phi_h, Phi_m, ratio, eca, eca0; 
//  double w, e;
public:
  double iT, m0, h0;
  double G_Ca;
  IT_TC(double v) {
     G_Ca = 2; //2;///////////////////////////////////
     Phi_m = pow(Qm,((Cels-24)/10));
     Phi_h = pow(Qh,((Cels-24)/10));
//     m0 = 1 / (1+exp(-(v+60.5)/6.2));
//     h0 = 1 / (1+exp((v+84)/4.03));
     m0 = 1 / (1+exp(-(v+59)/6.2));////////////////////////
     h0 = 1 / (1+exp((v+83)/4.0));
     eca0 = 1000*8.31441*(273.15 + Cels)/(2*96489);  } 
  void calc(double m, double h, double &fm, double &fh,  
            double v, double cai, double x);
};

double IT_TC::Shift = 2, IT_TC::Ca_0 = 2, IT_TC::Cels = 36;
double IT_TC::Qm = 3.55, IT_TC::Qh = 3; //2.8;

void IT_TC::calc(double m, double h, double &fm, double &fh,
                 double v, double cai, double x) {
  ratio = Ca_0/cai;
    if(ratio <= 0.)printf("\n LOG ERROR: RE: cai=%lf ratio=%lf",cai,ratio);
  eca = eca0 * log(ratio);
  iT = G_Ca*m*m*h*(v - eca); 

/*    w = v * 0.001 * 2 * 96480 / (8.314*(36 + 273.16));
      if (fabs(w) > 1e-4) e = w/(exp(w)-1);
      else e = 1 - w/2;
      iT = G_Ca*m*m*h* (-0.001) * 2 * 96480 * (Ca_0 - cai*exp(w)) * e; */
/*    m_inf = 1 / (1+exp(-(v+60.5)/6.2));
      h_inf = 1 / (1+exp((v+84)/4.03));*/

  m_inf = 1 / (1+exp(-(v+59)/6.2));//////////////////////////////////////
  h_inf = 1 / (1+exp((v+83)/4.)); /////////////////////////////////////

//  iT = G_Ca*m_inf*m_inf*h*(v - eca);

  tau_m = (1/(exp(-(v+131.6)/16.7)+exp((v+16.8)/18.2)) + 0.612) / Phi_m;
//  tau_m = 0.15*m_inf*(1.7+exp(-(v+30.8)/13.5));/////////////////////
  tau_h = (30.8 + (211.4 + exp((v + Shift + 113.2)/5))/
           (1+exp((v + Shift + 84)/3.2))) / Phi_h;

  /*  if (v<-80) tau_h = exp((v+467)/66.6) / Phi_h; 
      else tau_h = (exp(-(v+21.88)/10.52)+28) / Phi_h;  */  	  

  fm = -(1/tau_m)*(m - m_inf);                                
  fh = -(1/tau_h)*(h - h_inf);
}

//----------------- h-current (TC cell) -----------------------------------
class Ih_TC {
  static double E_h, cac, pc, Shift, Cels, k2, k4, nca, nexp, taum;
  double h_inf, tau_s, alpha, beta, k3p, cc, Phi;
public:
  double ih, p10, o10, o20;
  double G_h, k1ca, ginc;

  Ih_TC(double v, double cai) {
     G_h = 0.02; //////////////
     ginc = 1.5;  //1.5;///          To decrease after-burst depolarisation
     Phi = pow(3,((Cels-36)/10));
     h_inf = 1/(1 + exp((v + 75 - Shift)/5.5));
     tau_s = (taum + 1000 / (exp((v + 71.5 - Shift)/14.2) + 
                          exp(-(v + 89 - Shift)/11.6))) / Phi;
     alpha = h_inf/tau_s;
     beta = (1 - h_inf)/tau_s;
     p10 = 1/(1 + pow((cac/cai),nca));
     o10 = 1/(1 + beta/alpha + pow((p10/pc),nexp));
     o20 = pow((p10/pc),nexp) * o10;
  }
  void calc(double o1, double p1, double o2,  
                 double &fo1, double &fp1, double &fo2,  
                 double v, double cai, double x);
};

double Ih_TC::E_h = -40, Ih_TC::cac = 0.002, Ih_TC::pc = 0.01; 
double Ih_TC::Shift = 0, Ih_TC::Cels = 36;
double Ih_TC::k2 = 0.0004 /*0.0004*/; /////To decrease after-burst depolarisation
double Ih_TC::k4 = 0.001;
double Ih_TC::nca = 4, Ih_TC::nexp = 1, Ih_TC::taum = 20;

void Ih_TC::calc(double o1, double p1, double o2, double &fo1, double &fp1, 
                 double &fo2, double v, double cai, double x) {
  ih = G_h*(o1 + ginc * o2)*(v - E_h);
  h_inf = 1/(1 + exp((v + 75 - Shift)/5.5));
  tau_s = (taum + 1000 / (exp((v + 71.5 - Shift)/14.2) + 
                          exp(-(v + 89 - Shift)/11.6))) / Phi;
  alpha = h_inf/tau_s;
  beta = (1 - h_inf)/tau_s;
  k1ca = k2 * pow((cai/cac),nca);
  k3p = k4 * pow((p1/pc),nexp);
  fo1 = alpha * (1-o1-o2) - beta * o1 + k4 * o2 - k3p * o1;
  fp1 = k1ca * (1-p1) - k2 * p1 + k4 * o2 - k3p * o1;
  fo2 = k3p * o1 - k4 * o2;
}

//----------------------Potassium A-current (TC cell)------------------------
class IA_TC {
  static double E_K, Cels;     
  double m_inf, tau_m, h_inf, tau_h, Tad;                                 
public:
  double iA, m0, h0, G_A;
  IA_TC(double v) {
    G_A = 10; //2; //////////////////////////////////////////////
    Tad = pow(3,((Cels-23.5)/10));
    m0 = 1.0 / (1+exp(-(v+60)/8.5));
    h0 = 1.0/(1+exp((v+78)/6)); } 
  void calc(double m, double h, double &fm, double &fh, double v, double x);
};

double IA_TC::Cels = 36, IA_TC::E_K = -95;

void IA_TC::calc(double m, double h, double &fm, double &fh, double v, double x){
  iA = G_A*m*m*m*m*h*(v - E_K);
  tau_m = (1.0/( exp((v+35.82)/19.69)+exp(-(v+79.69)/12.7) ) +0.37) / Tad;
  m_inf = 1.0 / (1+exp(-(v+60)/8.5));
  tau_h = 1.0/((exp((v+46.05)/5)+exp(-(v+238.4)/37.45))) / Tad;
  if(v >= -63) 
      tau_h = 19.0/Tad;
  h_inf = 1.0/(1+exp((v+78)/6));
  fm = -(1/tau_m)*(m - m_inf);                                  
  fh = -(1/tau_h)*(h - h_inf);
}

//---------------------Hight-threshold Ca2+ current (CX cell)----------------
class IHVA_CX {
  static double Shift, Ca_0, Cels, Qm, Qh, E_Ca;
  double m_inf, tau_m, h_inf, tau_h, Phi_m, Phi_h, a, b, vm;
  double ratio, eca0, eca;
public:
  double iHVA, m0, h0;
  double G_HVA;
  IHVA_CX(double v) {
    G_HVA = 0.03;
    Phi_m = pow(Qm,((Cels-23)/10));
    Phi_h = pow(Qh,((Cels-23)/10));
    vm = v + Shift;
    a = 0.055*(-27 - vm)/(exp((-27-vm)/3.8) - 1);
    b = 0.94*exp((-75-vm)/17);
    m0 = a/(a+b);
    a = 0.000457*exp((-13-vm)/50);
    b = 0.0065/(exp((-vm-15)/28) + 1);
    h0 = a/(a+b);
//    eca0 = 1000*8.31441*(273.15 + Cels)/(2*96489);
    } 
  void calc(double m, double h, double &fm, double &fh, double v, double cai, double x);
};

double IHVA_CX::Shift = 0, IHVA_CX::Ca_0 = 2, IHVA_CX::E_Ca = 140; 
double IHVA_CX::Qm = 2.3, IHVA_CX::Qh = 2.3, IHVA_CX::Cels = 23; //36;

void IHVA_CX::calc(double m, double h, double &fm, double &fh,
                 double v, double cai, double x) {
//------------ECa is fixed (=140mV) instead of using Nerst eq.-------------
//  ratio = Ca_0/cai;
//  if(ratio <= 0.)printf("\n LOG ERROR: RE: cai=%lf ratio=%lf",cai,ratio);
//  eca = eca0 * log(ratio);

  iHVA = Phi_m * G_HVA * m*m*h * (v - E_Ca);
  vm = v + Shift;

  a = 0.055*(-37 - vm)/(exp((-37-vm)/3.8) - 1);
  b = 0.94*exp((-55-vm)/17);
  tau_m = (1/(a+b))/Phi_m;
  m_inf = a/(a+b);

  a = 0.000457*exp((-13-vm)/50);
  b = 0.0065/(exp((-vm-15)/28) + 1);
  tau_h = (1/(a+b))/Phi_h;
  h_inf = a/(a+b);

  fm = -(m - m_inf)/tau_m;                                  
  fh = -(h - h_inf)/tau_h;
}

//--------------Ca-dependent potassium current (CX cell)-----------------------
class IKCa_CX {
  static double E_KCa, Ra, Rb, Cels, Q, caix;     
  double m_inf, tau_m, Tad, a, b;                                 
public:
  double iKCa, m0;
  double G_KCa;
  IKCa_CX(double cai) {
    G_KCa = 0.3; 
    Tad = pow(Q,((Cels-23)/10));
    a = Ra * cai;  //------becouse caix = 1
    b = Rb;
    m0 = a/(a+b);
  }
  void calc(double m, double &fm, double v, double cai, double x);
};

double IKCa_CX::E_KCa = -90, IKCa_CX::Q = 2.3, IKCa_CX::caix = 1;
double IKCa_CX::Ra = 0.01, IKCa_CX::Rb = 0.02, IKCa_CX::Cels = 23; //36;   

void IKCa_CX::calc(double m, double &fm, double v, double cai, double x){
  iKCa = Tad * G_KCa * m * (v - E_KCa);                         

//  a = Ra * pow(cai,caix);
  a = Ra * cai;  //------becouse caix = 1
  b = Rb;
  tau_m = (1/(a+b))/Tad;
  m_inf = a/(a+b);
  fm = -(1/tau_m)*(m - m_inf);                   
}

//--------------------Potassium M-current (CX cell)-------------------------
class IKm_CX {
  static double E_Km, Ra, Rb, Cels, Q, tha, qa;     
  double m_inf, tau_m, Tad, a, b;                                 
public:
  double iKm, m0;
  double G_Km;
  IKm_CX(double v) {
    G_Km = 0.01; 
    Tad = pow(Q,((Cels-23)/10));
    a = Ra * (v - tha) / (1 - exp(-(v - tha)/qa));
    b = -Rb * (v - tha) / (1 - exp((v - tha)/qa));
    m0 = a/(a+b);
  }
  void calc(double m, double &fm, double v, double x);
};

double IKm_CX::E_Km = -90, IKm_CX::Q = 2.3;
double IKm_CX::tha = -30, IKm_CX::qa = 9;
double IKm_CX::Ra = 0.001, IKm_CX::Rb = 0.001, IKm_CX::Cels = 36;   

void IKm_CX::calc(double m, double &fm, double v, double x){
  iKm = Tad * G_Km * m * (v - E_Km);
  a = Ra * (v - tha) / (1 - exp(-(v - tha)/qa));
  b = -Rb * (v - tha) / (1 - exp((v - tha)/qa));
  tau_m = (1/(a+b))/Tad;
  m_inf = a/(a+b);
  fm = -(1/tau_m)*(m - m_inf);                   
}

//--------------------Fast potassium current (CX cell)-------------------
class IKv_CX {
  static double Ra, Rb, Cels, Q, tha, qa;     
  double m_inf, tau_m, Tad, a, b;                                 
public:
  static double E_Kv;
  double iKv, g_Kv, m0;
  double G_Kv;
  IKv_CX(double v) {
    G_Kv = 150; 
    Tad = pow(Q,((Cels-23)/10));

    a = Ra * (v - tha) / (1 - exp(-(v - tha)/qa));
    b = -Rb * (v - tha) / (1 - exp((v - tha)/qa));
    m0 = a/(a+b);
  }
  void calc(double m, double &fm, double v, double x);
};

double IKv_CX::E_Kv = -90, IKv_CX::Q = 2.3;
double IKv_CX::tha = 25, IKv_CX::qa = 9;
double IKv_CX::Ra = 0.02, IKv_CX::Rb = 0.002, IKv_CX::Cels = 36;   

void IKv_CX::calc(double m, double &fm, double v, double x){
  g_Kv = Tad * G_Kv * m;
  iKv = g_Kv * (v - E_Kv);                         
  a = Ra * (v - tha) / (1 - exp(-(v - tha)/qa));
  b = -Rb * (v - tha) / (1 - exp((v - tha)/qa));
  tau_m = (1/(a+b))/Tad;
  m_inf = a/(a+b);
  fm = -(1/tau_m)*(m - m_inf);                   
}

//----------------Fast sodium current (CX cell)--------------------------
class INa_CX {
  static double Shift, Ca_0, Cels, Qm, Qh;
  static double tha, qa, Ra, Rb, thi1, thi2, qi, thinf, qinf, Rg, Rd; 
  double m_inf, tau_m, h_inf, tau_h, Phi_m, Phi_h, a, b, vm;
  double trap0(double v, double th, double a, double q) {
        if (fabs(v/th) > 1.0e-6) {
                return ( a * (v - th) / (1 - exp(-(v - th)/q)) );
        } else {
	  return (a * q ); }
        }
public:
  static double E_Na;
  double iNa, g_Na, m0, h0;
  double G_Na;
  INa_CX(double v) {
    G_Na = 3000;
    Phi_m = pow(Qm,((Cels-23)/10));
    Phi_h = pow(Qh,((Cels-23)/10));
    vm = v + Shift;
    a = trap0(vm,tha,Ra,qa);
    b = trap0(-vm,-tha,Rb,qa);
    m0 = a/(a+b);

    a = trap0(vm,thi1,Rd,qi);
    b = trap0(-vm,-thi2,Rg,qi);
    h0 = 1/(1+exp((vm-thinf)/qinf));
    } 
  void calc(double m, double h, double &fm, double &fh, 
            double v, double x);
};

double INa_CX::Shift = -10, INa_CX::E_Na = 60; 
double INa_CX::Qm = 2.3, INa_CX::Qh = 2.3, INa_CX::Cels = 36;
double INa_CX::tha = -35, INa_CX::qa = 9;
double INa_CX::Ra = 0.182,INa_CX::Rb = 0.124; 
double INa_CX::thi1 = -50, INa_CX::thi2 = -75, INa_CX::qi = 5;
double INa_CX::thinf = -65, INa_CX::qinf = 6.2;
double INa_CX::Rg = 0.0091, INa_CX::Rd = 0.024; 

void INa_CX::calc(double m, double h, double &fm, double &fh,
                 double v, double x) {

  g_Na = Phi_m * G_Na * m*m*m*h;
  iNa = g_Na * (v - E_Na);
  vm = v + Shift;

  a = trap0(vm,tha,Ra,qa);
  b = trap0(-vm,-tha,Rb,qa);
  tau_m = (1/(a+b))/Phi_m;
  m_inf = a/(a+b);

                //"h" inactivation 
  a = trap0(vm,thi1,Rd,qi);
  b = trap0(-vm,-thi2,Rg,qi);
  tau_h = (1/(a+b))/Phi_h;
  h_inf = 1/(1+exp((vm-thinf)/qinf));

  fm = -(m - m_inf)/tau_m;                                  
  fh = -(h - h_inf)/tau_h;
}


//------------------Low-theshold Ca2+ current (LN cell)-----------------
class IT_LN {
  static double Ca_0, Cels, Qm, Qh, Shift, E_Ca;
  double m_inf, tau_m, h_inf, tau_h, Phi_h, Phi_m, ratio, eca, eca0; 
public:
  double iT, m0, h0;
  double G_Ca;
  IT_LN(double v) {
     G_Ca = 2; //2;///////////////////////////////////
     Phi_m = pow(Qm,((Cels-24)/10));
     Phi_h = pow(Qh,((Cels-24)/10));
     m0 = 1 / (1+exp(-(v+20)/6.5));////////////////////////
     h0 = 1 / (1+exp((v+25)/12.));
     eca0 = 1000*8.31441*(273.15 + Cels)/(2*96489);  } 
  void calc(double m, double h, double &fm, double &fh,  
            double v, double cai, double x);
};

double IT_LN::Shift = 2, IT_LN::Ca_0 = 2, IT_LN::Cels = 24; //36;
double IT_LN::Qm = 3.55, IT_LN::Qh = 3, IT_LN::E_Ca = 140;

void IT_LN::calc(double m, double h, double &fm, double &fh,
                 double v, double cai, double x) {
//  ratio = Ca_0/cai;
//  if(ratio <= 0.)printf("\n LOG ERROR: RE: cai=%lf ratio=%lf",cai,ratio);
//  eca = eca0 * log(ratio);
  eca = E_Ca;
  iT = G_Ca*m*m*h*(v - eca); 

  m_inf = 1 / (1+exp(-(v+20)/6.5));//////////////////////////////////////
  h_inf = 1 / (1+exp((v+25)/12)); /////////////////////////////////////

  tau_m = 1.5; //1 + (v+30)*0.014;
    //(1/(exp(-(v+131.6)/167.)+exp((v+16.8)/182.)) + 0.612) / Phi_m;
  tau_h = 0.3*exp((v-40)/13) + 0.002*exp(-(v-60)/29);

//               (30.8 + (211.4 + exp((v + Shift + 113.2)/5))/
//               (1+exp((v + Shift + 84)/3.2))) / Phi_h;

  fm = -(1/tau_m)*(m - m_inf);                                
  fh = -(1/tau_h)*(h - h_inf);
}

//--------------fast K current (LN cells)------------------
class IKv_LN {
  static double Cels;
  double Alpha3, Beta3, v2, Phi;
  double tau_n, n_inf;
public:
  static double E_K;
  double iKv, n0;
  double G_Kv, Vtr;
  double S2;
  IKv_LN(double v) {
    G_Kv = 10;///////////////////////
    Vtr = -50;
    S2 = 0.02;
    v2 = v - Vtr;
    Phi = pow(3,((Cels-36)/10));

    Alpha3 = 0.032*(15 - v2)/(exp((15 - v2)/5) - 1);
    Beta3 = 0.5*exp((10 - v2)/40);
    n0 = Alpha3/(Alpha3 + Beta3);     } 
  void calc(double n, double &fn, double v, double x);
};

double IKv_LN::E_K = -95, IKv_LN::Cels = 22; 

void IKv_LN::calc(double n, double &fn, double v, double x){

  v2 = v - Vtr;
  iKv = G_Kv* n*n*n*n*(v - E_K);    
  Alpha3 = S2*(15 - v2)/(exp((15 - v2)/5) - 1);
  Beta3 = 0.5*exp((10 - v2)/40);
  tau_n = 1/(Alpha3 + Beta3) / Phi;
  n_inf = Alpha3/(Alpha3 + Beta3);
  
  fn  = -(n - n_inf)/tau_n;                 
}

//===================Now we'll CREATE the CELLS==================================
//-------------------RE CELL-----------------------------------------------
class RE: public IT_RE, public IKCa_RE, public ICAN_RE, 
          public INaK, public ICa {
  static double Cai0, V0;
public:
  double G_kl, G_l, E_l, S;

  RE() :IT_RE(V0), IKCa_RE(Cai0), ICAN_RE(Cai0), INaK(V0), ICa() {
        E_l = -70;
        G_l = 0.05;
        G_kl = 0.018; //0.012; //0.015;
        S = 1.43e-4;
        }
  void init(double *y) {
        y[0] = V0;
	y[1] = Cai0;
        y[2] = IT_RE::m0;
        y[3] = IT_RE::h0;
//        y[4] = IKCa_RE::m0;
//        y[5] = ICAN_RE::m0;
        y[4] = INaK::m0;
        y[5] = INaK::h0;
        y[6] = INaK::n0;
        }
  void calc(double x, double *y, double *f); 
};  

double RE::V0 = -61, RE::Cai0 = 0.0001;

void RE::calc(double x, double *y, double *f){
    IT_RE::calc(y[2], y[3], f[2], f[3], y[0], y[1], x);
//    IKCa_RE::calc(y[4], f[4], y[0],y[1], x);
//    ICAN_RE::calc(y[5], f[5], y[0],y[1], x);
//    INaK::calc(y[6], y[7], y[8], f[6], f[7], f[8], y[0], x); 
    INaK::calc(y[4], y[5], y[6], f[4], f[5], f[6], y[0], x);
    ICa::calc(y[1], f[1], iT, x);
    f[0] = -G_l * (y[0] - E_l) - iT /*- iKCa - iCAN*/ - iNa - iK
                         - G_kl * (y[0] - INaK::E_K);
}

//-------------------TC CELL-----------------------------------------------
class TC: public IT_TC, public Ih_TC, public INaK, public ICa, public IA_TC {
  static double G_l, Cai0, V0;
public:
  double G_kl, S, E_l;

  TC() :Ih_TC(V0,Cai0), IT_TC(V0), INaK(V0), ICa(), IA_TC(V0) {
        G_kl = 0.012;
        E_l = -70;
        INaK::G_Na = 90;
        INaK::G_K = 10;
        S = 2.9e-4;
        }
  void init(double *y) {
        y[0] = V0;
	y[1] = Cai0;
        y[2] = IT_TC::m0;
        y[3] = IT_TC::h0;
        y[4] = Ih_TC::o10;
        y[5] = Ih_TC::p10;
        y[6] = Ih_TC::o20;
        y[7] = INaK::m0;
        y[8] = INaK::h0;
        y[9] = INaK::n0;
        y[10] = IA_TC::m0;
        y[11] = IA_TC::h0;
        }
  void calc(double x, double *y, double *f); 
};  

double TC::Cai0 = 0.0001;
double TC::G_l = 0.01;//0.05; //0.01;///
double TC::V0 = -68;

void TC::calc(double x, double *y, double *f){
    IT_TC::calc(y[2], y[3], f[2], f[3], y[0], y[1], x);
    Ih_TC::calc(y[4], y[5], y[6], f[4], f[5], f[6], y[0], y[1], x);
    INaK::calc(y[7], y[8], y[9], f[7], f[8], f[9], y[0], x); 
    IA_TC::calc(y[10], y[11], f[10], f[11], y[0], x);
    ICa::calc(y[1], f[1], iT, x);
    f[0] = -G_l * (y[0] - E_l) - iT - ih - iNa - iK - iA
              - G_kl * (y[0] - INaK::E_K);
}

//-------------------CX CELL------------------------------------------------
//-------------------CX CELL (DENDRITE)-------------------------------------
class CX_DEND: public IHVA_CX, public IKCa_CX, public IKm_CX, 
               public INa_CX, public ICa {
  static double E_l, G_l;
public:
  double iDEND, I_Stim1;
  CX_DEND(double V0, double Cai0) :IHVA_CX(V0), IKCa_CX(Cai0), IKm_CX(V0),
                                   INa_CX(V0), ICa() { 
        INa_CX::G_Na = 1.5; 
        I_Stim1 = 0; 
        ICa::Taur = 150; //200;
        }
  void init(double V0, double Cai0, double *y) {
        y[0] = V0;
	y[1] = Cai0;
        y[2] = IHVA_CX::m0;
        y[3] = IHVA_CX::h0;
        y[4] = IKCa_CX::m0;
        y[5] = IKm_CX::m0;
        y[6] = INa_CX::m0;
        y[7] = INa_CX::h0;
  }
  void calc(double x, double *y, double *f); 
};  
double CX_DEND::E_l = -70; double CX_DEND::G_l = 1.0e3/30000;  //  mS/cm^2

void CX_DEND::calc(double x, double *y, double *f){
    IHVA_CX::calc(y[2], y[3], f[2], f[3], y[0], y[1], x);
    IKCa_CX::calc(y[4], f[4], y[0], y[1], x);
    IKm_CX::calc(y[5], f[5], y[0], x);
    INa_CX::calc(y[6], y[7], f[6], f[7], y[0], x);
    ICa::calc(y[1], f[1], iHVA, x);
    iDEND =  -G_l * (y[0] - E_l) - iHVA - iKCa - iKm - iNa + I_Stim1;     
}

//-------------------CX CELL (SOMA)-------------------------------------
class CX_SOMA: public IKv_CX, public INa_CX {
public:
  double v_SOMA, iSOMA, g1_SOMA, g2_SOMA, I_Stim2;
  CX_SOMA(double V0, double Cai0) :IKv_CX(V0), INa_CX(V0) { 
        I_Stim2 = 0; }
  void init(double V0, double Cai0, double *y) {
        v_SOMA = V0;
        y[0] = IKv_CX::m0;
        y[1] = INa_CX::m0;
        y[2] = INa_CX::h0;
        }
  void calc(double x, double *y, double *f); 
};  
void CX_SOMA::calc(double x, double *y, double *f){
    IKv_CX::calc(y[0], f[0], v_SOMA, x);
    INa_CX::calc(y[1], y[2], f[1], f[2], v_SOMA, x);
    g1_SOMA = g_Na + g_Kv;
    g2_SOMA = g_Na * INa_CX::E_Na + g_Kv * IKv_CX::E_Kv + I_Stim2; 
    iSOMA =  - iNa - iKv;     
}

//------------CX CELL (connect DENDRITE and SOMA)---------------------------
class CX: public CX_DEND, public CX_SOMA {
  static double Cai0, V0, C;
public:
  double kappa, rho, S_CX_SOMA, S_CX_DEND; 
  CX() :CX_DEND(V0,Cai0), CX_SOMA(V0,Cai0) {
        kappa = 10.0e3;      // kOm: to get mS=1/kOm 
        rho = 165; //140; //200; //165;
        }
  void init(double *y) {
        CX_DEND::init(V0, Cai0, y);
	CX_SOMA::init(V0, Cai0, y+N_DEND);
        S_CX_SOMA = 1.0e-6;  // cm^2
        S_CX_DEND = S_CX_SOMA * rho;        
        }
  void calc(double x, double *y, double *f); 
};  

double CX::Cai0 = 0.0001, CX::V0 = -68;
double CX::C = 0.75;   // uF/cm^2

void CX::calc(double x, double *y, double *f){
    CX_SOMA::calc(x, y+N_DEND, f+N_DEND);
    v_SOMA = (y[0] + kappa * S_CX_SOMA * g2_SOMA) / 
                            (1 + kappa*S_CX_SOMA * g1_SOMA);
    CX_DEND::calc(x, y, f);
    f[0] = (1.0/C) * ( iDEND + 1.0/(kappa*S_CX_DEND) * (v_SOMA - y[0]) );
}

//-------------------INTERNEURON-----------------------------------------------
class IN: public INaK, public IA_TC  {
  static double V0;
public:
  double G_kl, G_l, E_l, S, DC;

  IN() :INaK(V0), IA_TC(V0) {
        E_l = -55;
        G_l = 0.15;
        G_kl = 0.05;
        DC = 0;
        INaK::G_Na = 50;
        INaK::G_K = 10;
        INaK::Vtr = -55;
        INaK::VtrK = -45;
        S2 = 0.03;
        S = 1.43e-4;
        }
  void init(double *y) {
        y[0] = V0;
        y[1] = INaK::m0;
        y[2] = INaK::h0;
        y[3] = INaK::n0;
        y[4] = IA_TC::m0;
        y[5] = IA_TC::h0;
        }
  void calc(double x, double *y, double *f); 
};  

double IN::V0 = -61;

void IN::calc(double x, double *y, double *f){
    INaK::calc(y[1], y[2], y[3], f[1], f[2], f[3], y[0], x);
    IA_TC::calc(y[4], y[5], f[4], f[5], y[0], x);
    f[0] = -G_l * (y[0] - E_l) - iNa - iK - iA -
           G_kl * (y[0] - INaK::E_K) + DC/S;
}


//-------------------LN cell-----------------------------------------------

class LN: public INaK, public IT_LN, public ICa, public IKCa_CX, 
          public IKv_LN {
  static double V0, Cai0;
public:
  double G_kl, G_l, E_l, S, DC;

  LN() :INaK(V0), IT_LN(V0), IKCa_CX(Cai0), ICa(), IKv_LN(V0) {
    E_l = -50; //-58; //-54;
        G_l = 0.15;
        G_kl = 0.02;
        DC = 0;
        ICa::Taur = 150;
        ICa::D = 0.25;
	//        IKv_LN::Vtr = -45;
        S = 1.43e-4;
        }
  void init(double *y) {
        y[0] = V0;	
        y[1] = Cai0;
        y[2] = IT_LN::m0;
        y[3] = IT_LN::h0;
        y[4] = IKCa_CX::m0;
        y[5] = IKv_LN::n0;
        y[6] = INaK::m0;
        y[7] = INaK::h0;
        y[8] = INaK::n0;
        }
  void calc(double x, double *y, double *f); 
};  

double LN::V0 = -61, LN::Cai0 = 0.0001;

void LN::calc(double x, double *y, double *f){
    INaK::calc(y[6], y[7], y[8], f[6], f[7], f[8], y[0], x);
    IT_LN::calc(y[2], y[3], f[2], f[3], y[0], y[1], x);
    IKCa_CX::calc(y[4], f[4], y[0], y[1], x);
    ICa::calc(y[1], f[1], iT, x);
    IKv_LN::calc(y[5], f[5], y[0], x);
    f[0] = -G_l * (y[0] - E_l)  - iNa - iK - iT - iKCa - iKv -
           G_kl * (y[0] - INaK::E_K) + DC/S;
}


//---------SYNAPCES DESCRIPTION-------------------------------------------
//---------first order kiner model for GABA-A synapce---------------------
class Gaba_A {
  static double Cdur, Cmax, Deadtime, Prethresh;  
  double R, C, R0, R1;
  double lastrelease;
  double q, Rinf, Rtau;
  double exptable(double z)
    {
    if((z > -10) && (z < 10)) return( exp(z) );
    else return( 0 );
    }
public:
  double I, E_GABA, Alpha, Beta;
  Gaba_A() {
    E_GABA = -70; //-75; (-70 is more realistic?)
    R = 0, C = 0, R0 = 0, R1 = 0;
    Alpha = 10.5;
    Beta = 0.166;
    lastrelease = -100;
    Rinf = Cmax*Alpha / (Cmax*Alpha + Beta);
    Rtau = 1 / ((Alpha * Cmax) + Beta);
  }
  void calc(double g_GABA_A, double x, double y_pre, double y_post);
}; 
double Gaba_A::Cdur = 0.3, Gaba_A::Cmax = 0.5, Gaba_A::Deadtime = 1;
double Gaba_A::Prethresh = -20;
void Gaba_A::calc(double g_GABA_A, double x, double y_pre, 
                    double y_post) {

        q = ((x - lastrelease) - Cdur);         
        if (q > Deadtime) {
                if (y_post > Prethresh) {        
                        C = Cmax;                
                        R0 = R;
                        lastrelease = 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)));
        }
       I = g_GABA_A * R * (y_pre - E_GABA);
}

//------second order kiner model (including G-proteins) for GABA-B synapce----
class Gaba_B {
  static double E_GABA, Cmax, Deadtime, Prethresh;   
  static double Kd, n; 
  double Gn, q;
public:
  double C, lastrelease;
  double I, r0, g0, Gn1, Cdur, K1, K2, K1K2, K3, K4, G;
  Gaba_B() {
    Cdur = 0.3; 
    K1 = 1.0; //0.52; 
//    K1K2 = 400;  //     K1K2 = K1/K2
    K2 = 0.0025; //0.0013;
    K3 = 0.1; //0.098;
    K4 = 0.06; //0.1; //0.033;
    lastrelease = -10000000;
    C = 0, r0 = 0, g0 = 0;
  }
  void calc(double r, double g, double &fr, double &fg, 
            double g_GABA_B, double x, double y_pre, double y_post);
};
double Gaba_B::E_GABA = -95, Gaba_B::Cmax = 0.5, Gaba_B::Deadtime = 1;
double Gaba_B::Prethresh = -20;
double Gaba_B::Kd = 100, Gaba_B::n = 4; 

//double Gaba_B::Prethresh = 0, Gaba_B::Kd = 100, Gaba_B::n = 4;
//double Gaba_B::K1 = 0.70; //0.63;//0.57;//0.8;//0.52;//0.09;//AAAAA 0.52; 
//double Gaba_B::K2 = 0.0035;//0.0025;//0.0022;//0.002;//0.0013;//0.0012;//AAAAA 0.0013;
//double Gaba_B::K3 = 0.15;   //0.098;//0.18;//AAAAAAAA 0.098; 
//double Gaba_B::K4 = 0.05;   //0.033;//0.034;//AAAAAAAA 0.033

void Gaba_B::calc(double r, double g, double &fr, double &fg, 
                  double g_GABA_B, double x, double y_pre, double y_post) {
        Gn = pow(g,n); 
        Gn1 = Gn/(Gn + Kd);
        G= g_GABA_B * Gn1;
        I = G * (y_pre - E_GABA);

        q = ((x - lastrelease) - Cdur);         
        if (q > Deadtime) {
                if (y_post > Prethresh) {        
                        C = Cmax;                
                        lastrelease = x; }
        } else if (q < 0) {                     
        } else if (C == Cmax) {   C = 0;  }
        fr = K1 * C * (1 - r) - r * K2;
        fg = K3 * r - K4 * g;
}

class GB: public Gaba_B {
public:
//  double y[N_GB], f[N_GB];
  GB() :Gaba_B() { }
  void init(double *y){
      lastrelease = -10000000;
      C = 0;
      y[0] = 0;
      y[1] = 0;
      }   
  void calc(double g_GABA_B, double x, double *y, double *f, 
                                             double y_pre, double y_post){
       Gaba_B::calc(y[0], y[1], f[0], f[1], g_GABA_B, x, y_pre, y_post); 
       } 
};  

//------second order kiner model (including G-proteins) for GABA-B synapce----
class Gaba_BF {
  static double E_GABA, Cmax, Deadtime, Prethresh;   
  static double Kd, n; 
  double Gn, q;
  double F, TrF, dF;
public:
  double C, lastrelease;
  double I, r0, g0, Gn1, Cdur, K1, K2, K1K2, K3, K4, G, FF;
  Gaba_BF() {
    Cdur = 0.3; 
    K1 = 1.0; //0.52; 
//    K1K2 = 400;  //     K1K2 = K1/K2
    K2 = 0.0025; //0.0013;
    K3 = 0.1; //0.098;
    K4 = 0.06; //0.1; //0.033;
    F=1;
    dF=0; //0.25;
    TrF=10000;
    lastrelease = -10000000;
    C = 0, r0 = 0, g0 = 0;
  }
  void calc(double r, double g, double &fr, double &fg, 
            double g_GABA_B, double x, double y_pre, double y_post);
};
double Gaba_BF::E_GABA = -95, Gaba_BF::Cmax = 0.5, Gaba_BF::Deadtime = 1;
double Gaba_BF::Prethresh = -20;
double Gaba_BF::Kd = 100, Gaba_BF::n = 4; 

//double Gaba_B::Prethresh = 0, Gaba_B::Kd = 100, Gaba_B::n = 4;
//double Gaba_B::K1 = 0.70; //0.63;//0.57;//0.8;//0.52;//0.09;//AAAAA 0.52; 
//double Gaba_B::K2 = 0.0035;//0.0025;//0.0022;//0.002;//0.0013;//0.0012;//AAAAA 0.0013;
//double Gaba_B::K3 = 0.15;   //0.098;//0.18;//AAAAAAAA 0.098; 
//double Gaba_B::K4 = 0.05;   //0.033;//0.034;//AAAAAAAA 0.033

void Gaba_BF::calc(double r, double g, double &fr, double &fg,  
                  double g_GABA_B, double x, double y_pre, double y_post) {
        Gn = pow(g,n); 
        Gn1 = Gn/(Gn + Kd);

        FF = g_GABA_B * F;
        G= g_GABA_B * Gn1;
        I = G * F * (y_pre - E_GABA);

        q = ((x - lastrelease) - Cdur);         
        if (q > Deadtime) {
                if (y_post > Prethresh) {        
                        C = Cmax;      
                        F = 1 + (F + dF - 1.0) * exp(-q/TrF);          
                        lastrelease = x; }
        } else if (q < 0) {                     
        } else if (C == Cmax) {   C = 0;  }
        fr = K1 * C * (1 - r) - r * K2;
        fg = K3 * r - K4 * g;
}

class GBF: public Gaba_BF {
public:
//  double y[N_GB], f[N_GB];
  GBF() :Gaba_BF() { }
  void init(double *y){
      lastrelease = -10000000;
      C = 0;
      y[0] = 0;
      y[1] = 0;
      }   
  void calc(double g_GABA_B, double x, double *y, double *f, 
                                             double y_pre, double y_post){
       Gaba_BF::calc(y[0], y[1], f[0], f[1], g_GABA_B, x, y_pre, y_post); 
       } 
};  

//------------first order kiner model for AMPA synapce---------------------
class AMPA {
  static double E_AMPA;
  static double Cdur, Cmax, Deadtime, Prethresh, Cdel; 
  static double Alpha, Beta; 
  int s;
  double R, C, R0, R1;
  double lastrelease, lastspike;
  double q, Rinf, Rtau;
  double exptable(double z)
    {
    if((z > -10) && (z < 10)) return( exp(z) );
    else return( 0 );
    }
public:
  double I;
  AMPA() {
    R = 0, C = 0, R0 = 0, R1 = 0;
    lastrelease = -100;
    lastspike = -100;
    s = 1;
    Rinf = Cmax*Alpha / (Cmax*Alpha + Beta);
    Rtau = 1 / ((Alpha * Cmax) + Beta);
  }
  void calc(double g_AMPA, double x, double y_pre, double y_post);
};
double AMPA::E_AMPA = 0, AMPA::Cdur = 0.3, AMPA::Cmax = 0.5, AMPA::Deadtime = 1;
double AMPA::Cdel = 6; //synaptic delay to get some latency between depolarization
                         //of the presynaptic membrain and transmitter reliase
double AMPA::Prethresh = 0, AMPA::Alpha = 1, AMPA::Beta = 0.2;
void AMPA::calc(double g_AMPA, double x, double y_pre, 
                    double y_post) {

        q = ((x - lastrelease) - Cdur); 
        
        if(q > Deadtime) {
               if(y_post > Prethresh) {
                  if( (x - lastspike) > (Cdel + Cdur) ){
                     lastspike = x;
                     s = 1; } }  //the flag that spike was but wasn't utilized yet

               if((s == 1) && ((x - lastspike) > Cdel)) {
                  s = 0;         //spike was utilized
                  C = Cmax;                
                  R0 = R;
                  lastrelease = 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)));
        }
       I = g_AMPA * R * (y_pre - E_AMPA);
}

//------------first order kiner model for AMPA synapce with Facilitation-------
class AMPA_F {
  static double E_AMPA;
  static double Cdur, Cmax, Deadtime, Prethresh, Cdel; 
  static double Alpha, Beta; 
  int s;
  double R, C, R0, R1;
  double lastrelease, lastspike;
  double q, Rinf, Rtau;
  double F, dF, TrF;
  double exptable(double z)
    {
    if((z > -10) && (z < 10)) return( exp(z) );
    else return( 0 );
    }
public:
  double I, FF;
  AMPA_F() {
    R = 0, C = 0, R0 = 0, R1 = 0;
    lastrelease = -100;
    lastspike = -100;
    s = 1;
    F=1;
    dF=0; //0.25;
    TrF=10000;
    Rinf = Cmax*Alpha / (Cmax*Alpha + Beta);
    Rtau = 1 / ((Alpha * Cmax) + Beta);
  }
  void calc(double g_AMPA, double x, double y_pre, double y_post);
};
double AMPA_F::E_AMPA = 0, AMPA_F::Cdur = 0.3, AMPA_F::Cmax = 0.5, AMPA_F::Deadtime = 1;
double AMPA_F::Cdel = 6; //synaptic delay to get some latency between depolarization
                         //of the presynaptic membrain and transmitter reliase
double AMPA_F::Prethresh = 0, AMPA_F::Alpha = 1, AMPA_F::Beta = 0.2;

void AMPA_F::calc(double g_AMPA, double x, double y_pre, 
                    double y_post) {

        q = ((x - lastrelease) - Cdur); 
        
        if(q > Deadtime) {
               if(y_post > Prethresh) {
                  if( (x - lastspike) > (Cdel + Cdur) ){
                     lastspike = x;
                     s = 1; } }  //the flag that spike was but wasn't utilized yet

               if((s == 1) && ((x - lastspike) > Cdel)) {
                  s = 0;         //spike was utilized
                  C = Cmax;                
                  R0 = R;
                  F = 1; // + (F + dF - 1.0) * exptable(-q/TrF);
                  lastrelease = 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)));
        }
       FF = g_AMPA * F;
       I = FF * R * (y_pre - E_AMPA);
}


//------------first order kiner model for AMPA synapce WITH depression--------------
class AMPA_D1 {
  static double E_AMPA;
  static double Cdur, Cmax, Deadtime, Prethresh, Cdel; 
  static double Alpha, Beta; 
  int s;
  double R, C, R0, R1;
  double lastrelease, lastspike, E, Use, Tr;
  double q, Rinf, Rtau;
  double exptable(double z)
    {
    if((z > -10) && (z < 10)) return( exp(z) );
    else return( 0 );
    }
public:
  double I;
  AMPA_D1() {
    R = 0, C = 0, R0 = 0, R1 = 0;
    lastrelease = -100;
    lastspike = -100;
    s = 1;
    Rinf = Cmax*Alpha / (Cmax*Alpha + Beta);
    Rtau = 1 / ((Alpha * Cmax) + Beta);
    E = 1;
    Use = 0.5; 
    Tr = 500;
  }
  void calc(double g_AMPA, double x, double y_pre, double y_post);
};
double AMPA_D1::E_AMPA = 0, AMPA_D1::Cdur = 0.3, AMPA_D1::Cmax = 0.5, AMPA_D1::Deadtime = 1;
double AMPA_D1::Cdel = 0; //synaptic delay to get some latency between depolarization
                         //of the presynaptic membrain and transmitter reliase
double AMPA_D1::Prethresh = 0, AMPA_D1::Alpha = 0.94, AMPA_D1::Beta = 0.18;
void AMPA_D1::calc(double g_AMPA, double x, double y_pre, 
                    double y_post) {

        q = ((x - lastrelease) - Cdur); 
        
        if(q > Deadtime) {
               if(y_post > Prethresh) {
                  if( (x - lastspike) > (Cdel + Cdur) ){
                     lastspike = x;
                     s = 1; } }  //the flag that spike was but wasn't utilized yet

               if((s == 1) && ((x - lastspike) > Cdel)) {
                  s = 0;         //spike was utilized
                  C = Cmax;                
                  R0 = R;
                  E = 1 - (1 - E*(1-Use)) * exptable(-q/Tr);
                  lastrelease = 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)));
        }
       I = g_AMPA * R * E * (y_pre - E_AMPA);
}

//---------first order kiner model for GABA-A synapce with DEPRESSION---------------------
class Gaba_A_D1 {
  static double Cdur, Cmax, Deadtime, Prethresh; 
  static double Alpha, Beta;  
  double R, C, R0, R1;
  double lastrelease, E, Use, Tr;
  double q, Rinf, Rtau;
  double exptable(double z)
    {
    if((z > -10) && (z < 10)) return( exp(z) );
    else return( 0 );
    }
public:
  double I, E_GABA;
  Gaba_A_D1() {
    E_GABA = -70; //-75; (-70 is more realistic?)
    R = 0, C = 0, R0 = 0, R1 = 0;
    lastrelease = -100;
    Rinf = Cmax*Alpha / (Cmax*Alpha + Beta);
    Rtau = 1 / ((Alpha * Cmax) + Beta);
    E = 1;
    Use = 0.5; 
    Tr = 500;
  }
  void calc(double g_GABA_A, double x, double y_pre, double y_post);
}; 
double Gaba_A_D1::Cdur = 0.3, Gaba_A_D1::Cmax = 0.5, Gaba_A_D1::Deadtime = 1;
double Gaba_A_D1::Prethresh = 0, Gaba_A_D1::Alpha = 10.5, Gaba_A_D1::Beta = 0.166; 
void Gaba_A_D1::calc(double g_GABA_A, double x, double y_pre, 
                    double y_post) {

        q = ((x - lastrelease) - Cdur);         
        if (q > Deadtime) {
                if (y_post > Prethresh) {        
                        C = Cmax;                
                        R0 = R;
                        E = 1 - (1 - E*(1-Use)) * exptable(-q/Tr);
                        lastrelease = 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)));
        }
       I = g_GABA_A * E * R * (y_pre - E_GABA);
}

//-----first order kiner model for AMPA synapce 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 = 10, w=0.01, wom=0;
  }
  void init(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;
}

//-----first order kiner model for AMPA synapce used for external stimulation----
class Extern_ampa1 {
  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, sig, Noise;
  Extern_ampa1() {
    R = 0, C = 0, R0 = 0, R1 = 0, sig = 0.8;
    Noise = 0.86;
  }
  void init(unsigned int seek) {srand(seek);}
  void calc(double g_Extern_ampa, double x);
};
void Extern_ampa1::calc(double g_Extern_ampa, double x) {
       RRR = 2.0*rand()/(RAND_MAX + 1.0) - 1.0;
       R = 1.0 + RRR/sig + 0.5*(Noise-0.86);
       g = g_Extern_ampa * R;
}

//-----first order kiner model for external pulse----
class Pulse {
  static double Cdur, Cmax, Deadtime; 
  double R, C, R0, R1, Rinf1, T1;
  double lastrelease;
  double q, Rinf, Rtau;
  double TR;
  double exptable(double z)
    {
    if((z > -10) && (z < 10)) return( exp(z) );
    else return( 0 );
    }
public:
  double g, Alpha, Beta, Beta1, RS;
  Pulse() {
    Alpha = 0.01;
    Beta = 0.005;
    Beta1 = 1./120.; //1./150.;
    T1 = 250; //400;
    R = 0, C = 0, R0 = 0, R1 = 0;
    lastrelease = -1500;
    TR = 2000.0;
    Rinf = Cmax*Alpha / (Cmax*Alpha + Beta);
    Rtau = 1 / ((Alpha * Cmax) + Beta);
    Rinf1 = 1/(1+ exp(Beta1 * (-T1)));
  }
  void init() {
    Rinf = Cmax*Alpha / (Cmax*Alpha + Beta);
    Rtau = 1 / ((Alpha * Cmax) + Beta);
    }
  void calc(double x);
};
double Pulse::Cdur = 500, Pulse::Cmax = 1, Pulse::Deadtime = 1;
void Pulse::calc(double x) {

        q = ((x - lastrelease) - Cdur);         
        if (q > Deadtime) {
                if ( (x - lastrelease) > TR ) {        
                        C = Cmax;                
                        R0 = R;
                        lastrelease = 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)));
//           R = (R1/Rinf1)/(1+ exp (Beta1 * (x-(lastrelease+Cdur)-T1)));
        }
     RS=R/Rinf;
}

//---------first order kiner model for gradient GABA-A synapce--------------
class Gaba_Ag {
  static double Cdur, Cmax, Deadtime, Prethresh, Sig;  
  double C;
  double exptable(double z)
    {
    if((z > -10) && (z < 10)) return( exp(z) );
    else return( 0 );
    }
public:
  double I, E_GABA, Alpha, Beta, G;
  Gaba_Ag() {
    E_GABA = -70; //-75; (-70 is more realistic?)
    C = 0;
    Alpha = 10;
    Beta = 0.2;
  }
  void calc(double g_GABA_A, double x, double *r, double *fr,
                               double y_pre, double y_post);
}; 
//double Gaba_Ag::Cdur = 0.3, Gaba_Ag::Cmax = 0.5, Gaba_Ag::Deadtime = 1;
double Gaba_Ag::Prethresh = -20, Gaba_Ag::Sig = 1.5;
void Gaba_Ag::calc(double g_GABA_A, double x, double *y, double *f, 
                    double y_post, double y_pre) {
   
        G = g_GABA_A * y[0];
        I = G * (y_post - E_GABA);

        C = 1/(1+exp(-(y_pre - Prethresh)/Sig));
        f[0] = Alpha * C * (1 - y[0]) - y[0] * Beta;

}

//---------first order kiner model for gradient GABA-A synapce FACILIT--------
class Gaba_AgF {
  static double Cdur, Cmax, Deadtime, Prethresh, Sig;  
  double C, F, dF, TrF, q;
  double exptable(double z)
    {
    if((z > -10) && (z < 10)) return( exp(z) );
    else return( 0 );
    }
public:
  double I, E_GABA, Alpha, Beta, FF, lastrelease, G;
  Gaba_AgF() {
    E_GABA = -70; //-75; (-70 is more realistic?)
    C = 0;
    F=1;
    dF=0; //0.2; //0.25;
    TrF=10000;
    lastrelease = -10000000;
    Alpha = 10;
    Beta = 0.2;
  }
  void calc(double g_GABA_A, double x, double *r, double *fr,
                               double y_pre, double y_post);
}; 
double Gaba_AgF::Cdur = 10, Gaba_AgF::Cmax = 0.5, Gaba_AgF::Deadtime = 1;
double Gaba_AgF::Prethresh = -20, Gaba_AgF::Sig = 1.5;
void Gaba_AgF::calc(double g_GABA_A, double x, double *y, double *f, 
                    double y_post, double y_pre) {

        FF = g_GABA_A * F;
        G = FF * y[0];
        I = G * (y_post - E_GABA);

        C = 1/(1+exp(-(y_pre - Prethresh)/Sig));
        q = ((x - lastrelease) - Cdur);         
        if (q > Deadtime) {
                if (y_pre > Prethresh) {        
                        F = 1 + (F + dF - 1.0) * exptable(-q/TrF);
                        lastrelease = x;
                        }
                } 
        f[0] = Alpha * C * (1 - y[0]) - y[0] * Beta;
}


//------second order kiner model (including G-proteins) for GABA-B synapce---
class Gaba_Bg {
  static double E_GABA, Cmax, Deadtime, Prethresh, Sig;   
  static double Kd, n; 
  double Gn, q;
public:
  double C, lastrelease;
  double I, r0, g0, Gn1, Cdur, K1, K2, K1K2, K3, K4;
  Gaba_Bg() {
    Cdur = 0.3; 
    K1 = 0.52; 
    K2 = 0.0013; 
    K3 = 0.098;
    K4 = 0.033;
    lastrelease = -10000000;
    C = 0, r0 = 0, g0 = 0;
  }
  void calc(double r, double g, double &fr, double &fg, 
            double g_GABA_B, double x, double y_pre, double y_post);
};
double Gaba_Bg::E_GABA = -95, Gaba_Bg::Cmax = 0.5, Gaba_Bg::Deadtime = 1;
double Gaba_Bg::Prethresh = -20,  Gaba_Bg::Sig = 0.1;
double Gaba_Bg::Kd = 100, Gaba_Bg::n = 4; 

void Gaba_Bg::calc(double r, double g, double &fr, double &fg, 
                  double g_GABA_B, double x, double y_pre, double y_post) {
        Gn = pow(g,n); 
        Gn1 = Gn/(Gn + Kd);
        I = g_GABA_B * Gn1 * (y_pre - E_GABA);

        C = 1/(1+exp(-(y_post - Prethresh)/Sig));
       
        fr = K1 * C * (1 - r) - r * K2;
        fg = K3 * r - K4 * g;
}

class GBg: public Gaba_Bg {
public:
//  double y[N_GB], f[N_GB];
  GBg() :Gaba_Bg() { }
  void init(double *y){
      lastrelease = -10000000;
      C = 0;
      y[0] = 0;
      y[1] = 0;
      }   
  void calc(double g_GABA_B, double x, double *y, double *f, 
                                             double y_pre, double y_post){
       Gaba_Bg::calc(y[0], y[1], f[0], f[1], g_GABA_B, x, y_pre, y_post); 
       } 
};  

//+++++++++++++++++++ MAIN PROGRAM +++++++++++++++++++++++++++++++++++++++++++

//----------external functuons----------------------------------------------
void rk(unsigned, void (unsigned, double, double*, double*), double, double, 
                                double*, double*, double*, double*);
void fun(unsigned, double, double*, double*);
void solout (long, double, double, double*, unsigned, int*);
//----------external variables---------------------------------------------
double g_gaba_a, g_gaba_a1, g_gaba_b, g_gaba_b_re_re, g_ampa, g_ampa1;
double g_ampa_cx_cx,  g_ampa_cx_tc, g_ampa_cx_re, g_ampa_tc_cx; 
double g_ampa_cx_in,  g_gaba_a_in_cx, g_gaba_b_in_cx, g_ampa_tc_in; 
double g_ext_tc, g_ext_re, g_ext_cx, g_ext_in;
double g_ext1, g_ext2, g_ext3, g_ext4, Noise, tt;
double SynapsesTCRE[Mre*Mre1*N_TC_RE][5], SynapsesTCTC[Mtc*Mtc1*N_TC_TC][5], SynapsesRETC[Mtc*Mtc1*N_RE_TC][5],  SynapsesRERE[Mre*Mre1*N_RE_RE][5], SynapsesRETCB[Mtc*Mtc1*N_RE_TC][5]; 

FILE *f1, *f2, *f3, *f4, *f6, *f7, *f8, *f9, *f10, *f11, *f12, *f13, *f22, *f222, *fNoise, *fSYNre, *fSYNtc, *fsRERE, *fsRETC, *fsTCRE, *fsTCTC, *fsRETCB, *fre, *ftc,*fC_PN,*fC_LN,*fC_LNLN,*fC_LNPN,*fC_PNLN,*fC_PNPN, *ftest1, *ftest2, *fStimPN, *fStimLN,*fGABA_A_TCIP, *fGABA_B_TCIP,*fAMPA_TCIP;

double xx,xx1,xx2,xx3,xx4,xx5;

int no_re[Mre][Mre1][N_RE], no_tc[Mtc][Mtc1][N_TC];
int no_g[Mtc][Mtc1][N_RE_TC][N_GB], no_gb_rere[Mre][Mre1][N_RE_RE][N_GB];
int no_cx[Mc][Mc1][N_CX], no_in[Mc][Mc1][N_IN], no_gcx[Mc][Mc1][N_IN_CX][N_GB];
int no_ga_retc[Mtc][Mtc1][N_RE_TC][N_GA], no_ga_rere[Mre][Mre1][N_RE_RE][N_GA];
double C_RERE[Mre][Mre1][Mre][Mre1], C_RETC[Mre][Mre1][Mtc][Mtc1];
double C_TCRE[Mtc][Mtc1][Mre][Mre1], C_TCTC[Mtc][Mtc1][Mtc][Mtc1];
double C_RE[Mre][Mre1], C_TC[Mtc][Mtc1], C_RE1[Mre][Mre1], C_TC1[Mtc][Mtc1];
int k_re_re[Mre][Mre1], k_re_tc[Mtc][Mtc1], k_tc_re[Mre][Mre1], k_tc_tc[Mtc][Mtc1];
double inh_a[Mtc][Mtc1],inh_b[Mtc][Mtc1],exc_AMPA[Mtc][Mtc1];
//int **C_RE, **C_TC;
int KMAXrere, KMAXretc, KMAXtcre, KMAXtctc, KMAXretcB;
/*
//----------external variables ---------------------------------------------
double g_gaba_a, g_gaba_a1, g_gaba_b, g_gaba_b_re_re, g_ampa, g_ampa1;
double g_ampa_cx_cx,  g_ampa_cx_tc, g_ampa_cx_re, g_ampa_tc_cx; 
double g_ampa_cx_in,  g_gaba_a_in_cx, g_gaba_b_in_cx, g_ampa_tc_in; 
double g_ext_tc, g_ext_re, g_ext_cx, g_ext_in;
double g_ext1, g_ext2, g_ext3, g_ext4, Noise, tt;
double SynapsesTCRE[Mre*Mre1*N_TC_RE][5], SynapsesTCTC[Mtc*Mtc1*N_TC_TC][5], SynapsesRETC[Mtc*Mtc1*N_RE_TC][5],  SynapsesRERE[Mre*Mre1*N_RE_RE][5], SynapsesRETCB[Mtc*Mtc1*N_RE_TC][5]; 

FILE *f1, *f2, *f3, *f4, *f6, *f7, *f8, *f9, *f10, *f11, *f12, *f13, *f22, *f222, *fNoise, *fSYNre, *fSYNtc, *fsRERE, *fsRETC, *fsTCRE, *fsTCTC, *fsRETCB, *fre, *ftc;

int no_re[Mre][Mre1][N_RE], no_tc[Mtc][Mtc1][N_TC];
int no_g[Mtc][Mtc1][N_RE_TC][N_GB], no_gb_rere[Mre][Mre1][N_RE_RE][N_GB];
int no_cx[Mc][Mc1][N_CX], no_in[Mc][Mc1][N_IN], no_gcx[Mc][Mc1][N_IN_CX][N_GB];
int no_ga_retc[Mtc][Mtc1][N_RE_TC][N_GA], no_ga_rere[Mre][Mre1][N_RE_RE][N_GA];
long int C_RERE[Mre][Mre1][Mre][Mre1], C_RETC[Mre][Mre1][Mtc][Mtc1];
long int C_TCRE[Mtc][Mtc1][Mre][Mre1], C_TCTC[Mtc][Mtc1][Mtc][Mtc1];
long int C_RE[Mre][Mre1], C_TC[Mtc][Mtc1], C_RE1[Mre][Mre1], C_TC1[Mtc][Mtc1];
long int k_re_re[Mre][Mre1], k_re_tc[Mtc][Mtc1], k_tc_re[Mre][Mre1], k_tc_tc[Mtc][Mtc1];
//int **C_RE, **C_TC;
long int KMAXrere, KMAXretc, KMAXtcre, KMAXtctc, KMAXretcB;
*/
//----------external classes (beginning of initialization)------------------
Gaba_AgF      *g_re_re[Mre][Mre1];
Gaba_AgF      *g_re_tc[Mtc][Mtc1];
GBF          *gb[Mtc][Mtc1];
GBF          *gb_re_re[Mre][Mre1];
AMPA_F         *a_tc_re[Mre][Mre1];
AMPA_F         *a_tc_tc[Mtc][Mtc1];
LN           re_cell[Mre][Mre1];
IN           tc_cell[Mtc][Mtc1];

AMPA         *a_cx_cx[Mc][Mc1];
AMPA         *a_cx_in[Mc][Mc1];
Gaba_A       *ga_in_cx[Mc][Mc1];
GB           *gb_in_cx[Mc][Mc1];
CX           cx_cell[Mc][Mc1];
CX           in_cell[Mc][Mc1];         

AMPA         *a_cx_tc[Mtc][Mtc1];
AMPA         *a_cx_re[Mre][Mre1];
AMPA         *a_tc_cx[Mc][Mc1];
AMPA         *a_tc_in[Mc][Mc1];

Extern_ampa1  a_ext1, a_ext2, a_ext3, a_ext4;
Pulse aP, aPRE[Mre][Mre1];

//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
int main(int argc,char **argv)
{
//---------allocate place for ALL variables and functions-----------
  double y_ini[N_EQ], f_ini[N_EQ];

//---------allocate place for TWO temporal arrays using by RK.C solver 
  double y1[N_EQ], y2[N_EQ];

//---------general parameters----------------------------------------------
  printf("\n ... reading input parameters");
  double t = 0, tmax, tmax1, t3D, ttime, TAU, *Noi;
  double h = 0.04, red, t_ext, g_Extern_ampa3, g_Extern_ampa4, R, r[2];
  int i, j, k, l, i1, j1, ii = 0, i_inp, i_gip = 0, ih, ni;
  double g_GABA_A, g_GABA_A1, g_GABA_B, g_AMPA, g_AMPA1, g_GABA_B_RE_RE;
  double g_AMPA_CX_CX, g_AMPA_CX_TC, g_AMPA_CX_RE, g_AMPA_TC_CX;
  double g_AMPA_CX_IN, g_GABA_A_IN_CX, g_GABA_B_IN_CX, g_AMPA_TC_IN;
  double g_Extern_ampa1, g_Extern_ampa2, g_Extern_ampa, av=0, avr=0;
  int k1=0, p=0, pmax=0, nt=17501;
  int *vm[Mtc][Mtc1];
  double tsp[Mtc][Mtc1], *v[Mtc][Mtc1], *vs[Mtc][Mtc1];
  unsigned ic=10;
  int tc;
//---------solver parameters (DOPRI)-------------------------------------------
  int      res, iout = 2, itoler = 0; 
  double   atoler = 1.0E-5, rtoler = 1.0E-5; 

//----------arrays initialization----------------------------------------------
  for(i=0; i<N_EQ; i++){
    y_ini[i] = 0, f_ini[i] = 0; 
    y1[i] = 0, y2[i] = 0; }

//-------parameter initialization (from file)----------------------------------
if (argc <= 1) {
  puts("Command parameters");
  puts("-----------------------");
  puts("Input File"); }

if (!(f1=fopen(argv[1],"r"))) {
   printf("%s doesn't exist\n",argv[1]);
   exit(0); }

  fscanf(f1, "%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %d", 
      &tmax, &t3D, &ttime, &t_ext, 
      &g_GABA_A, &g_GABA_A1, &g_GABA_B, &g_GABA_B_RE_RE, &g_AMPA, &g_AMPA1,
      &g_AMPA_CX_CX, &g_AMPA_CX_TC, &g_AMPA_CX_RE, &g_AMPA_TC_CX, 
      &g_AMPA_CX_IN, &g_GABA_A_IN_CX, &g_GABA_B_IN_CX, &g_AMPA_TC_IN,  
      &g_Extern_ampa1, &g_Extern_ampa2, &g_Extern_ampa3, &g_Extern_ampa4, 
      &i_inp);
  printf("\n param: %lf %lf %lf %lf A=%lf A1=%lf B=%lf B1=%lf AM=%lf AM1=%lf AM_CX_CX=%lf AM_CX_TC=%lf AM_CX_RE=%lf AM_TC_CX=%lf AM_CX_IN=%lf GA_IN_CX=%lf GB_IN_CX=%lf AM_TC_IN=%lf  %lf %lf %lf %lf %d", 
       tmax, t3D, ttime, t_ext, 
       g_GABA_A, g_GABA_A1, g_GABA_B, g_GABA_B_RE_RE, g_AMPA, g_AMPA1, 
       g_AMPA_CX_CX, g_AMPA_CX_TC, g_AMPA_CX_RE, g_AMPA_TC_CX,
       g_AMPA_CX_IN, g_GABA_A_IN_CX, g_GABA_B_IN_CX, g_AMPA_TC_IN,  
       g_Extern_ampa1, g_Extern_ampa2, g_Extern_ampa3, g_Extern_ampa4, 
       i_inp);
  fclose(f1);

  cout << "\n RE-RE_MAX=" << MS_RE_RE_MAX << " RE-TC_MAX=" << MS_RE_TC_MAX << 
                                          " TC-RE_MAX=" << MS_TC_RE_MAX;

//----------classes initialization (continue)----------------------------
  for(i=0; i < Mre; ++i)
  for(j=0; j < Mre1; ++j){
    g_re_re[i][j] = new Gaba_AgF[N_RE_RE];
    a_tc_re[i][j] = new AMPA_F[N_TC_RE];
    a_cx_re[i][j] = new AMPA[N_CX_RE]; 
    if(I_GB_RERE == 1) gb_re_re[i][j] = new GBF[N_RE_RE];
  }

  for(i=0; i < Mtc; ++i)
  for(j=0; j < Mtc1; ++j){
    g_re_tc[i][j] = new Gaba_AgF[N_RE_TC];
    if(I_GB == 1) gb[i][j] = new GBF[N_RE_TC];
    a_tc_tc[i][j] = new AMPA_F[N_TC_TC];
    a_cx_tc[i][j] = new AMPA[N_CX_TC];
  }

  for(i=0; i < Mc; ++i)
  for(j=0; j < Mc1; ++j){ 
    if(I_CX == 1) a_cx_cx[i][j] = new AMPA[N_CX_CX];
    if(I_CX == 1 && I_IN == 1) a_cx_in[i][j] = new AMPA[N_CX_IN];
    if(I_CX == 1 && I_IN == 1) ga_in_cx[i][j] = new Gaba_A[N_IN_CX];
    if(I_CX == 1 && I_IN == 1) gb_in_cx[i][j] = new GB[N_IN_CX];

    if(I_CX == 1 && I_TC == 1) a_tc_cx[i][j] = new AMPA[N_TC_CX];
    if(I_IN == 1 && I_TC == 1) a_tc_in[i][j] = new AMPA[N_TC_IN];
  }

   for(i=0; i < Mc; ++i)
   for(j=0; j < Mc1; ++j)
     in_cell[i][j].rho = 50;         //

   for(i=0; i < Mtc; ++i) 
      for(j=0; j < Mtc1; ++j){
            v[i][j] = new double[nt];
            vs[i][j] = new double[nt];
            vm[i][j] = new int[nt];
}
   for(i=0; i < Mtc; ++i) 
      for(j=0; j < Mtc1; ++j)
         for(k=0; k < nt; ++k){
            v[i][j][k] = 0.0;
            vm[i][j][k] = 0;
            vs[i][j][k] = 0.0;
         }

   Noi = new double[2500000];

/* for(i=0; i < Mre; ++i) 
      for(j=0; j < Mre1; ++j){
         for(k=0; k < Mre; ++k) 
            C_RERE[i][j][k] = new int[Mre1];
         for(k=0; k < Mtc; ++k) 
            C_RETC[i][j][k] = new int[Mtc1];
}
   for(i=0; i < Mtc; ++i) 
      for(j=0; j < Mtc1; ++j){
         for(k=0; k < Mre; ++k) 
            C_TCRE[i][j][k] = new int[Mre1];
         for(k=0; k < Mtc; ++k) 
            C_TCTC[i][j][k] = new int[Mtc1];
}
   C_TC = (int**) new int[Mtc];
   for(i=0; i < Mtc; ++i) 
      C_TC[i] = new int[Mtc1];
   for(i=0; i < Mtc; ++i) 
      for(j=0; j < Mtc1; ++j)
          printf("\n TC %d %d %d", i,j,C_TC[i][j]);

   C_RE = (int**) new int[Mre];
   for(i=0; i < Mre; ++i) 
      C_RE[i] = new int[Mre1];
   for(i=0; i < Mre; ++i) 
      for(j=0; j < Mre1; ++j)
          printf("\n RE %d %d %d", i,j,C_RE[i][j]);
*/

//----------creating the integer arrays containing the addresses--------
//----------of  ALL internal variables for ALL objects RE, TC ----------
//----------and GB classes (e.g., no_re[i][j][k] is the address --------
//----------of the variable y[k] for the object re_cell[i][j]) ---------
//----------NOTE: this is the relative adresses and you should ---------
//----------add the REAL address of the first element of the -----------
//----------original 1D array to use these arrays-----------------------
  for(i=0; i < Mre; ++i)
  for(j=0; j < Mre1; ++j)
  for(k=0; k < N_RE; ++k)
     no_re[i][j][k] = k + (j + i*Mre1) * N_RE;
  for(i=0; i < Mtc; ++i)
  for(j=0; j < Mtc1; ++j)
  for(k=0; k < N_TC; ++k)
     no_tc[i][j][k] = Mre*Mre1*N_RE*I_RE + k + (j + i*Mtc1) * N_TC;  
  for(i=0; i < Mtc; ++i)
  for(j=0; j < Mtc1; ++j)
  for(k=0; k < N_RE_TC; ++k)
  for(l=0; l < N_GB; ++l)
     no_g[i][j][k][l] = Mre*Mre1*N_RE*I_RE + Mtc*Mtc1*N_TC*I_TC + l + 
                                   (k + (j + i*Mtc1)*N_RE_TC) * N_GB;
  for(i=0; i < Mtc; ++i)
  for(j=0; j < Mtc1; ++j)
  for(k=0; k < N_RE_TC; ++k)
  for(l=0; l < N_GA; ++l)
     no_ga_retc[i][j][k][l] = Mre*Mre1*N_RE*I_RE + Mtc*Mtc1*N_TC*I_TC + 
                              Mtc*Mtc1*N_RE_TC*N_GB*I_RE*I_TC*I_GB +    
                              l + (k + (j + i*Mtc1)*N_RE_TC) * N_GA;
  for(i=0; i < Mre; ++i)
  for(j=0; j < Mre1; ++j)
  for(k=0; k < N_RE_RE; ++k)
  for(l=0; l < N_GA; ++l)
     no_ga_rere[i][j][k][l] = Mre*Mre1*N_RE*I_RE + Mtc*Mtc1*N_TC*I_TC + 
                              Mtc*Mtc1*N_RE_TC*N_GB*I_RE*I_TC*I_GB + 
                              Mtc*Mtc1*N_RE_TC*N_GA*I_RE*I_TC +
                              l + (k + (j + i*Mre1)*N_RE_RE) * N_GA;
  for(i=0; i < Mre; ++i)
  for(j=0; j < Mre1; ++j)
  for(k=0; k < N_RE_RE; ++k)
  for(l=0; l < N_GB; ++l)
     no_gb_rere[i][j][k][l] = Mre*Mre1*N_RE*I_RE + Mtc*Mtc1*N_TC*I_TC +
                              Mtc*Mtc1*N_RE_TC*N_GB*I_RE*I_TC*I_GB + 
                              Mtc*Mtc1*N_RE_TC*N_GA*I_RE*I_TC +
                              Mre*Mre1*N_RE_RE*N_GA*I_RE +
                              + l + (k + (j + i*Mre1)*N_RE_RE) * N_GB;

  for(i=0; i < Mc; ++i)
  for(j=0; j < Mc1; ++j)
  for(k=0; k < N_CX; ++k)
     no_cx[i][j][k] = Mre*Mre1*N_RE*I_RE + Mtc*Mtc1*N_TC*I_TC + 
                              Mtc*Mtc1*N_RE_TC*N_GB*I_RE*I_TC*I_GB +
                              Mtc*Mtc1*N_RE_TC*N_GA*I_RE*I_TC + 
                              Mre*Mre1*N_RE_RE*N_GA*I_RE +
                              Mre*Mre1*N_RE_RE*N_GB*I_RE*I_GB_RERE +
                                    k + (j + i*Mc1) * N_CX;
  for(i=0; i < Mc; ++i)
  for(j=0; j < Mc1; ++j)
  for(k=0; k < N_IN; ++k)
     no_in[i][j][k] = Mre*Mre1*N_RE*I_RE + Mtc*Mtc1*N_TC*I_TC + 
                               Mtc*Mtc1*N_RE_TC*N_GB*I_RE*I_TC*I_GB +
                               Mtc*Mtc1*N_RE_TC*N_GA*I_RE*I_TC + 
                               Mre*Mre1*N_RE_RE*N_GA*I_RE +
                               Mre*Mre1*N_RE_RE*N_GB*I_RE*I_GB_RERE +
                               Mc*Mc1*N_CX*I_CX + k + (j + i*Mc1) * N_IN;
  for(i=0; i < Mc; ++i)
  for(j=0; j < Mc1; ++j)
  for(k=0; k < N_IN_CX; ++k)
  for(l=0; l < N_GB; ++l)
     no_gcx[i][j][k][l] = Mre*Mre1*N_RE*I_RE + Mtc*Mtc1*N_TC*I_TC + 
                          Mtc*Mtc1*N_RE_TC*N_GB*I_RE*I_TC*I_GB +
                          Mtc*Mtc1*N_RE_TC*N_GA*I_RE*I_TC + 
                          Mre*Mre1*N_RE_RE*N_GA*I_RE +
                          Mre*Mre1*N_RE_RE*N_GB*I_RE*I_GB_RERE +
                          Mc*Mc1*N_CX*I_CX + Mc*Mc1*N_IN*I_IN +
                          l + (k + (j + i*Mc1)*N_IN_CX) * N_GB;
                                   
//---variable initialization (additional for standart constructor)----------
  for(i=0; i<Mre; i++)
    for(j=0; j<Mre1; j++){
      if(I_RE == 1) re_cell[i][j].init(y_ini+no_re[i][j][0]);
    }

  for(i=0; i<Mtc; i++)
    for(j=0; j<Mtc1; j++){
      if(I_TC == 1) tc_cell[i][j].init(y_ini+no_tc[i][j][0]);
//      if(I_RE == 1 && I_TC == 1 && I_GB == 1) for(k=0; k < N_RE_TC; k++){
//                           gb[i][j][k].init(y_ini+no_g[i][j][k][0]); }    
    }

  for(i=0; i<Mc; i++)
    for(j=0; j<Mc1; j++){
      if(I_CX == 1) cx_cell[i][j].init(y_ini+no_cx[i][j][0]); 
      if(I_IN == 1) in_cell[i][j].init(y_ini+no_in[i][j][0]);
      if(I_CX == 1 && I_IN == 1) for(k=0; k < N_IN_CX; k++){
                       gb_in_cx[i][j][k].init(y_ini+no_gcx[i][j][k][0]); }     
    }
  a_ext1.init(8);
  a_ext2.init(9);
  a_ext3.init(10);
  a_ext4.init(11);
    for(j = 0; j < Mtc1; ++j)
       for(i = 0; i < Mtc; ++i) 
          tsp[i][j] = -100; 


printf("\n NOISE");
fNoise = fopen("time_Hz100Syn200", "r");
i=0;
    while(tt< 99000.){
       fscanf(fNoise,"%lf %lf",&tt, &Noise);
       Noi[i] = Noise;
       //       printf("\n %lf %lf",tt,Noise);
       i=i+1;
    }
fclose(fNoise);
//----------here we are changing some variables----------------------


  for(i=0; i < Mre; ++i)
  for(j=0; j < Mre1; ++j){
     re_cell[i][j].G_kl = 0.01; //0.0015; //0.005;
  }
  for(i=0; i < Mtc; ++i)
  for(j=0; j < Mtc1; ++j){
    tc_cell[i][j].G_kl = 0.04; //0.01; //0.012; //0.015;
  }
  for(i=0; i < Mre; ++i)
  for(j=0; j < Mre1; ++j){
     re_cell[i][j].G_KCa = 0.25;
     //     re_cell[i][j].G_Ca = 0;
     re_cell[i][j].G_Kv = 7;
  }
  a_ext1.sig = 1000000000.0;
  a_ext2.sig = 1000000000.0;

  //  a_ext3.sig = 1000000000.0;
  //  a_ext4.sig = 1000000000.0;
//---changes of the parameters to get variability------------------------
     r[0] = 2.0 * rand()/(RAND_MAX + 1.0) - 1.0;
     r[1] = 2.0 * rand()/(RAND_MAX + 1.0) - 1.0;

srand(1);
  if(I_RE == 1)
  for(i=0; i < Mre; ++i)
   for(j=0; j < Mre1; ++j){
A:   R = 2.0 * rand()/(RAND_MAX + 1.0) - 1.0;
     if((R < -1) || (R > 1)) cout << "\n CHECK RANDOM !!!!!!!!";
//     if(R < -0.6) goto A;
//     if((r[0]*r[1]>0) && (r[1]*R>0)) goto A;
     while((r[0]*r[1]>0) && (r[1]*R>0)) 
                         {R = 2.0 * rand()/(RAND_MAX + 1.0) - 1.0;}
     r[1]=r[0]; r[0]=R;
     re_cell[i][j].G_kl = re_cell[i][j].G_kl + R * 0.005; //0.015; //0.02; 
     cout << "\n R=" << R << "  G_l(RE)=" << re_cell[i][j].G_kl; }
cout << "\n -------------------------------------------------";

srand(3);
  if(I_TC == 1) 
  for(i=0; i < Mtc; ++i)
   for(j=0; j < Mtc1; ++j){
     R = 2.0 * rand()/(RAND_MAX + 1.0) - 1.0;
     if((R < -1) || (R > 1)) cout << "\n CHECK RANDOM !!!!!!!!";
     while((r[0]*r[1]>0) && (r[1]*R>0)) {R = 2.0 * rand()/(RAND_MAX + 1.0) - 1.0;}
     r[1]=r[0]; r[0]=R;
     tc_cell[i][j].G_kl = tc_cell[i][j].G_kl + R * 0.005; 
     cout << "\n R=" << R << "  G_kl(TC)=" << tc_cell[i][j].G_kl; }
cout << "\n -------------------------------------------------";

/* srand(2);
   if(I_TC == 1)
   for(i=0; i < Mtc; ++i)
    for(j=0; j < Mtc1; ++j){
     R = 2.0 * rand()/(RAND_MAX + 1.0) - 1.0;
     if((R < -1) || (R > 1)) cout << "\n CHECK RANDOM !!!!!!!!";
     while((r[0]*r[1]>0) && (r[1]*R>0)) {R = 2.0 * rand()/(RAND_MAX + 1.0) - 1.0;}
     r[1]=r[0]; r[0]=R;
     tc_cell[i][j].G_h = tc_cell[i][j].G_h + R * 0.002; 
     cout << "\n R=" << R << "  G_h(TC)=" << tc_cell[i][j].G_h; }
cout << "\n -------------------------------------------------";
*/

for(i=0; i<Mtc; i++)
for(j=0; j<Mtc1; j++){
  R = 1.0 * rand()/(RAND_MAX + 1.0);
  if(R < 0.1) { tc_cell[i][j].DC=0.0001;
                printf("\n TC %d %d", i,j);}
  }

for(i=0; i<Mre; i++)
for(j=0; j<Mre1; j++){
  R = 1.0 * rand()/(RAND_MAX + 1.0);
  if(R < 0.1) { re_cell[i][j].DC=0.0001;
                printf("\n RE %d %d", i,j);}
  }

//--------------end variability------------------------------------
printf("\n Zero connections");
srand(7);
for(i=0; i<Mtc; i++)
for(j=0; j<Mtc1; j++) {
   printf("\n %d %d %d", i,j,C_TC[i][j]);
   C_TC[i][j]=0;
}

for(i=0; i<Mre; i++)
for(j=0; j<Mre1; j++) {
   printf("\n %d %d %d", i,j,C_RE[i][j]);
   C_RE[i][j]=0;
}

for(i=0; i<Mre; i++)
for(j=0; j<Mre1; j++)
for(i1=0; i1<Mre; i1++)
for(j1=0; j1<Mre1; j1++) C_RERE[i][j][i1][j1]=0;

for(i=0; i<Mre; i++)
for(j=0; j<Mre1; j++)
for(i1=0; i1<Mtc; i1++)
for(j1=0; j1<Mtc1; j1++) C_RETC[i][j][i1][j1]=0;

for(i=0; i<Mtc; i++)
for(j=0; j<Mtc1; j++)
for(i1=0; i1<Mre; i1++)
for(j1=0; j1<Mre1; j1++) C_TCRE[i][j][i1][j1]=0;

for(i=0; i<Mtc; i++)
for(j=0; j<Mtc1; j++)
for(i1=0; i1<Mtc; i1++)
for(j1=0; j1<Mtc1; j1++) C_TCTC[i][j][i1][j1]=0;

printf("\n Start connections \n");

fC_PN = fopen("C_PN.txt","r");

for(i=0; i<Mtc; i++)
for(j=0; j<Mtc1; j++){
   fscanf(fC_PN,"%lf",&xx);
   C_TC[i][j] = xx;
// printf("\n C_TC=%lf",C_TC[i][j]);
    printf("%lf\n",C_TC[i][j]);
  }
 fclose(fC_PN);

fC_LN = fopen("C_LN.txt","r");

for(i=0; i<Mre; i++)
for(j=0; j<Mre1; j++){
   fscanf(fC_LN,"%lf",&xx5);
   C_RE[i][j] = xx5;
   printf("%lf\n",C_RE[i][j]);
  }
 fclose(fC_LN);


fC_LNPN = fopen("C_LNPN.txt","r");

 for(i1 = 0; i1 < Mtc; i1++) 
   for(j1 = 0; j1 < Mtc1; j1++)
     for( i = 0; i < Mre; i++)
       for( j = 0; j < Mre1; j++){
	 fscanf(fC_LNPN,"%lf",&xx1);
	 C_RETC[i][j][i1][j1] = xx1;
//  printf("\n C_RETC=%lf",C_RETC[i][j][i1][j1]);
       }

// for(i1=0; i1<Mre; i1++)
// for(j1=0; j1<Mre1; j1++)
// for(i=0; i<Mtc; i++)
// for(j=0; j<Mtc1; j++){
//    fscanf(fC_LNPN,"%lf",&xx1);
//    C_RETC[i][j][i1][j1] = xx1;
// //    printf("\n CRETC=%lf",C_RETC[i][j][i1][j1]);
//   }
 fclose(fC_LNPN);

fC_LNLN = fopen("C_LNLN.txt","r");

for(i1 = 0; i1 < Mre; i1++)
  for(j1 = 0; j1 < Mre1; j1++)
    for(i=0; i<Mre; i++)
      for(j=0; j<Mre1; j++){
	fscanf(fC_LNLN,"%lf",&xx2);
	C_RERE[i][j][i1][j1] = xx2;
//  printf("\n C_RERE=%lf",C_RERE[i][j][i1][j1]);
      }
 fclose(fC_LNLN);


fC_PNLN = fopen("C_PNLN.txt","r");

for(i1=0; i1<Mre; i1++)
  for(j1=0; j1<Mre1; j1++)
    for(i=0; i<Mtc; i++)
      for(j=0; j<Mtc1; j++){
	fscanf(fC_PNLN,"%lf",&xx3);
	C_TCRE[i][j][i1][j1] = xx3;
// printf("\n C_TCRE=%lf",C_TCRE[i][j][i1][j1]);
      }
 fclose(fC_PNLN);

fC_PNPN = fopen("C_PNPN.txt","r");

for(i1=0; i1<Mtc; i1++)
  for(j1=0; j1<Mtc1; j1++)
    for(i=0; i<Mtc; i++)
      for(j=0; j<Mtc1; j++){
	fscanf(fC_PNPN,"%lf",&xx4);
	C_TCTC[i][j][i1][j1] = xx4;
 	printf("\n C_TCTC=%lf",C_TCTC[i][j][i1][j1]);
      }
 fclose(fC_PNPN);

/*
for(i=0; i<Mtc; i++)
for(j=0; j<Mtc1; j++){
  R = 1.0 * rand()/(RAND_MAX + 1.0);
  if(R < 0.3) {C_TC[i][j]=1; 
                printf("TC: %d %d \n", i,j);}
  }

for(i=0; i<Mre; i++)
for(j=0; j<Mre1; j++){
  R = 1.0 * rand()/(RAND_MAX + 1.0);
  if(R < 0.3) {C_RE[i][j]=1; 
                printf("RE: %d %d \n", i,j);}
  }

  
for(i=0; i<Mre; i++)
for(j=0; j<Mre1; j++)
for(i1=0; i1<Mre; i1++)
for(j1=0; j1<Mre1; j1++){
  R = 1.0 * rand()/(RAND_MAX + 1.0);
  if(R < 0.5) {C_RERE[i][j][i1][j1]=1;
//               C_RERE[i1][j1][i][j]=1;
}
  }

for(i=0; i<Mre; i++)
for(j=0; j<Mre1; j++)
for(i1=0; i1<Mtc; i1++)
for(j1=0; j1<Mtc1; j1++){
  R = 1.0 * rand()/(RAND_MAX + 1.0);
  if(R < 0.5) {C_RETC[i][j][i1][j1]=1;
//               C_TCRE[i1][j1][i][j]=1;
}
  }

for(i=0; i<Mtc; i++)
for(j=0; j<Mtc1; j++)
for(i1=0; i1<Mre; i1++)
for(j1=0; j1<Mre1; j1++){
  R = 1.0 * rand()/(RAND_MAX + 1.0);
  if(R < 0.5) {C_TCRE[i][j][i1][j1]=1;
//               C_RETC[i1][j1][i][j]=1;
}
  }

for(i=0; i<Mtc; i++)
for(j=0; j<Mtc1; j++)
for(i1=0; i1<Mtc; i1++)
for(j1=0; j1<Mtc1; j1++){
  R = 1.0 * rand()/(RAND_MAX + 1.0);
  if(R < 0.5) {C_TCTC[i][j][i1][j1]=1;
//               C_TCTC[i1][j1][i][j]=1;
}
  }
*/
/*
C_RETC[0][0][0][0]=1;
C_RETC[0][0][1][0]=1;
C_RETC[0][0][2][0]=1;
C_RETC[1][0][3][0]=1;
C_RETC[1][0][4][0]=1;
C_RETC[1][0][5][0]=1;

C_RERE[0][0][1][0]=1;
C_RERE[1][0][0][0]=1; 

C_RE[0][0]=1;
C_RE[1][0]=1;
C_TC[1][0]=1;
C_TC[3][0]=1;
*/

//--- Change Input ----------------------------------
/* srand(12);
 for(i=0; i<Mre; i++)
 for(j=0; j<Mre1; j++) {
    printf("\n %d %d %d", i,j,C_RE[i][j]);
    C_RE[i][j]=0;
 }
 for(i=0; i<Mre; i++)
 for(j=0; j<Mre1; j++){
   R = 1.0 * rand()/(RAND_MAX + 1.0);
   if(R < 0.3) {C_RE[i][j]=1;
                 printf("RE: %d %d \n", i,j);}
   }
*/

printf("\n End connections");
j=Mtc1/2;
for(i=0; i<Mtc; i++)
  printf("\n TC %ld %ld %ld", i, j, C_TC[i][j]);
j=Mtc1/4;
for(i=0; i<Mtc; i++)
  printf("\n TC %ld %ld %ld", i, j, C_TC[i][j]);
j=3*Mtc1/4;
for(i=0; i<Mtc; i++)
  printf("\n TC %ld %ld %ld", i, j, C_TC[i][j]);

for(i=0; i<Mre; i++)
  printf("\n RE %ld %ld", i, C_RE[i][Mre1/2]);

for(i1=0; i1<Mtc; i1++){
  for(i=0; i<Mre; i++) printf("%d ", C_RETC[i][0][i1][0]);
  printf("\n");
}

//------------------stimulus initialization----------------------
aP.init();
aPRE[0][0].Beta=0.0025; //0.0015;
srand(2);
   for(i=0; i < Mre; ++i)
    for(j=0; j < Mre1; ++j){
     R = 2.0 * rand()/(RAND_MAX + 1.0) - 1.0;
     if((R < -1) || (R > 1)) cout << "\n CHECK RANDOM !!!!!!!!";
     while((r[0]*r[1]>0) && (r[1]*R>0)) {R = 2.0 * rand()/(RAND_MAX + 1.0) - 1.0;}
     r[1]=r[0]; r[0]=R;
     aPRE[i][j].Beta = aPRE[0][0].Beta + R * 0.0001; 
     aPRE[0][0].Beta = 0.0025; //0.0015;
     aPRE[i][j].init();
}

//--------------open ALL files-------------------------------------
  f2 = fopen("dat", "w");
  f22 = fopen("dat1", "w");
  f222 = fopen("dat1all", "w");
  f6 = fopen("graf_re", "w");
  f7 = fopen("time_re", "w");
  f8 = fopen("graf_tc", "w");
  f9 = fopen("time_tc", "w");
  f10 = fopen("graf_cx", "w");
  f11 = fopen("time_cx", "w");
  f12 = fopen("graf_in", "w");
  f13 = fopen("time_in", "w");
  fNoise = fopen("time_Hz100Syn200", "r");
  fre = fopen("LNvariations","w");
  ftc = fopen("PNvariations","w");

//---------------show initial values of ALL variables----------------
  for(j = 0; j < Mtc1; ++j)
  for(i = 0; i < Mtc; ++i){
     if(I_TC == 1){
        printf("\n TC: ");
        for(k = 0; k < N_TC; ++k) 
           printf("%lf ", y_ini[no_tc[i][j][k]]);}
  }
  for(j = 0; j < Mre1; ++j)
  for(i = 0; i < Mre; ++i){
     if(I_RE == 1){
        printf("\n RE: ");
        for(k = 0; k < N_RE; ++k)   
           printf("%lf ", y_ini[no_re[i][j][k]]);}
  }
  for(j = 0; j < Mc1; ++j)
  for(i = 0; i < Mc; ++i){
     if(I_CX == 1){
        printf("\n CX: ");
        for(k = 0; k < N_CX; ++k) 
           printf("%lf ", y_ini[no_cx[i][j][k]]);}
     if(I_IN == 1){
        printf("\n IN: ");
        for(k = 0; k < N_IN; ++k)   
           printf("%lf ", y_ini[no_in[i][j][k]]);}
  }

  fSYNre = fopen("SYNre", "w"); 
  fSYNtc = fopen("SYNtc", "w");
  
  for(i1=0; i1<Mtc; i1++){
    fprintf(fSYNre,"%d ", i1);
    for(i=0; i<Mre; i++)
      fprintf(fSYNre,"%d ", C_RETC[i][0][i1][0]);
    fprintf(fSYNre,"\n");
  }
  for(i1=0; i1<Mtc; i1++){
    fprintf(fSYNtc,"%d ", i1);
    for(i=0; i<Mtc; i++)
      fprintf(fSYNtc,"%d ", C_TCTC[i][0][i1][0]);
    fprintf(fSYNtc,"\n");
  }
  
  fclose(fSYNre);
  fclose(fSYNtc);

//------------------------------------

 for(i=0; i<Mtc; i++)
   for(j=0; j<Mtc1; j++) 
     C_TC1[i][j]=C_TC[i][j];
 
 for(i=0; i<Mre; i++)
   for(j=0; j<Mre1; j++) 
     C_RE1[i][j]=C_RE[i][j];

tc=100000000;

//----------------CALCULATION----------------------------------------
  printf("\n CALCULATION IN PROGRESS!!!: t= %lf: tmax= %lf", t,tmax);
  ih = (int) (1/h);
  p=0;

  while( t < tmax){ 
  ii = ii + 1;       

  rk(N_EQ, fun, h, t, y_ini, f_ini, y1, y2);
  t = t + h;

if((t > 49.0)&&(t > 50.0)){
//--------the scaling of the conductances-----------------------------------
  //  g_ext3 = g_Extern_ampa/100;
  //  g_ext2 = g_Extern_ampa/100;
  g_ext1 = g_Extern_ampa1;
  g_ext2 = g_Extern_ampa2; 
  g_ext3 = g_Extern_ampa3;
  g_ext4 = g_Extern_ampa4;

if(I_RE == 1){
  g_gaba_a = g_GABA_A / re_cell[0][0].S; 
  g_ampa = g_AMPA / re_cell[0][0].S;
  g_gaba_b_re_re = g_GABA_B_RE_RE / re_cell[0][0].S;
//  g_ampa_cx_re = g_AMPA_CX_RE / re_cell[0][0].S;
  g_ext_re = g_Extern_ampa/5; ///re_cell[0][0].S;
}
if(I_TC == 1){
  g_gaba_a1 = g_GABA_A1 / tc_cell[0][0].S;
  g_ampa1 = g_AMPA1 / tc_cell[0][0].S;
  g_gaba_b = g_GABA_B / tc_cell[0][0].S;
//  g_ampa_cx_tc = g_AMPA_CX_TC / tc_cell[0][0].S;
  g_ext_tc = g_Extern_ampa; ///tc_cell[0][0].S;
}

if(I_CX == 1){
  g_ampa_cx_cx = g_AMPA_CX_CX / cx_cell[0][0].S_CX_DEND;
  g_ampa_tc_cx = g_AMPA_TC_CX / cx_cell[0][0].S_CX_DEND;
  g_gaba_a_in_cx = g_GABA_A_IN_CX / cx_cell[0][0].S_CX_DEND;
  g_gaba_b_in_cx = g_GABA_B_IN_CX / cx_cell[0][0].S_CX_DEND;
  g_ext_cx = g_Extern_ampa/cx_cell[0][0].S_CX_DEND/10; //2;
}
if(I_IN == 1){
  g_ampa_cx_in = g_AMPA_CX_IN / in_cell[0][0].S_CX_DEND;  //S_IN;
  g_ampa_tc_in = g_AMPA_TC_IN / in_cell[0][0].S_CX_DEND;  //S_IN;
//  g_ext_in = g_Extern_ampa/in_cell[0][0].S_CX_DEND / 3; // S_IN /15; //3;
}
}

Noise=Noi[ii];
//--- Change Input ----------------------------------

if((t > (tc-h/2.)) && (t < (tc+h/2.))){
 tc=tc+5000000;
 ic=ic+1;
 srand(ic);

 
 for(i=0; i<Mtc; i++)
   for(j=0; j<Mtc1; j++) {
     C_TC[i][j]=C_TC1[i][j];
 }

 for(i=0; i<Mre; i++)
   for(j=0; j<Mre1; j++) {
     C_RE[i][j]=C_RE1[i][j];
 }
/*
 fprintf(ftc,"t=%lf \n", t);
 for(i=0; i<Mtc; i++)
 for(j=0; j<Mtc1; j++){
   R = 1.0 * rand()/(RAND_MAX + 1.0);
   if((C_TC1[i][j]==0) && (R < 0.04)){
      C_TC[i][j]=1;
      fprintf(ftc,"TC+ %d %d \n", i,j);
   }}
 for(i=0; i<Mtc; i++)
 for(j=0; j<Mtc1; j++){
   R = 1.0 * rand()/(RAND_MAX + 1.0);
   if((C_TC1[i][j]==1) && (R < 0.1)){
      C_TC[i][j]=0;
      fprintf(ftc,"TC- %d %d \n", i,j);
   }}
*/
 fprintf(fre,"t=%lf \n", t);
 for(i=0; i<Mre; i++)
 for(j=0; j<Mre1; j++){
   R = 1.0 * rand()/(RAND_MAX + 1.0);
   if((C_RE1[i][j]==0) && (R < 0.2)){
      C_RE[i][j]=1;
      fprintf(fre,"RE+ %d %d \n", i);
   }}
/*
 for(i=0; i<Mre; i++)
 for(j=0; j<Mre1; j++){
   R = 1.0 * rand()/(RAND_MAX + 1.0);
   if((C_RE1[i][j]==1) && (R < 0.1)){
      C_RE[i][j]=0;
      fprintf(fre,"RE- %d %d \n", i);
   }}
*/
}


//-------------------------------------------------
if( (t>(1200-0.5*h)) && (t<(1200+0.5*h)) ){
KMAXrere=0;
KMAXtcre=0;
KMAXretc=0;
KMAXtctc=0;
KMAXretcB=0;
for(j = 0; j < Mre1; ++j)
  for(i = 0; i < Mre; ++i){
    for(k = 0; k < k_re_re[i][j]; ++k){
       KMAXrere=KMAXrere+1;
       SynapsesRERE[KMAXrere-1][0]=g_re_re[i][j][k].FF*re_cell[0][0].S;    //SynapsesRERE[KMAXrere-1][0]=g_re_re[i][j][k].FF;
    }
    for(k = 0; k < k_tc_re[i][j]; ++k){
       KMAXtcre=KMAXtcre+1;
       SynapsesTCRE[KMAXtcre-1][0]=a_tc_re[i][j][k].FF;
    }
  }
for(j = 0; j < Mtc1; ++j)
  for(i = 0; i < Mtc; ++i){
    for(k = 0; k < k_re_tc[i][j]; ++k){
       KMAXretc=KMAXretc+1;
       SynapsesRETC[KMAXretc-1][0]=g_re_tc[i][j][k].FF*tc_cell[0][0].S;         //SynapsesRETC[KMAXretc-1][0]=g_re_tc[i][j][k].FF;
    }
    for(k = 0; k < k_re_tc[i][j]; ++k){
       KMAXretcB=KMAXretcB+1;
       SynapsesRETCB[KMAXretcB-1][0]=gb[i][j][k].FF;
    }
    for(k = 0; k < k_tc_tc[i][j]; ++k){
       KMAXtctc=KMAXtctc+1;
       SynapsesTCTC[KMAXtctc-1][0]=a_tc_tc[i][j][k].FF;
    }
  }
}

if( (t>(3200-0.5*h)) && (t<(3200+0.5*h)) ){
KMAXrere=0;
KMAXtcre=0;
KMAXretc=0;
KMAXtctc=0;
KMAXretcB=0;
for(j = 0; j < Mre1; ++j)
  for(i = 0; i < Mre; ++i){
    for(k = 0; k < k_re_re[i][j]; ++k){
       KMAXrere=KMAXrere+1;
       SynapsesRERE[KMAXrere-1][1]=g_re_re[i][j][k].FF*re_cell[0][0].S;      //SynapsesRERE[KMAXrere-1][1]=g_re_re[i][j][k].FF;
    }
    for(k = 0; k < k_tc_re[i][j]; ++k){
       KMAXtcre=KMAXtcre+1;
       SynapsesTCRE[KMAXtcre-1][1]=a_tc_re[i][j][k].FF;
    }
  }
for(j = 0; j < Mtc1; ++j)
  for(i = 0; i < Mtc; ++i){
    for(k = 0; k < k_re_tc[i][j]; ++k){
       KMAXretc=KMAXretc+1;
       SynapsesRETC[KMAXretc-1][1]=g_re_tc[i][j][k].FF*tc_cell[0][0].S;        //SynapsesRETC[KMAXretc-1][1]=g_re_tc[i][j][k].FF;
    }
    for(k = 0; k < k_re_tc[i][j]; ++k){
       KMAXretcB=KMAXretcB+1;
       SynapsesRETCB[KMAXretcB-1][1]=gb[i][j][k].FF;
    }
    for(k = 0; k < k_tc_tc[i][j]; ++k){
       KMAXtctc=KMAXtctc+1;
       SynapsesTCTC[KMAXtctc-1][1]=a_tc_tc[i][j][k].FF;
    }
  }
}

if( (t>(5200-0.5*h)) && (t<(5200+0.5*h)) ){
KMAXrere=0;
KMAXtcre=0;
KMAXretc=0;
KMAXtctc=0;
KMAXretcB=0;
for(j = 0; j < Mre1; ++j)
  for(i = 0; i < Mre; ++i){
    for(k = 0; k < k_re_re[i][j]; ++k){
       KMAXrere=KMAXrere+1;
       SynapsesRERE[KMAXrere-1][2]=g_re_re[i][j][k].FF*re_cell[0][0].S;      //SynapsesRERE[KMAXrere-1][1]=g_re_re[i][j][k].FF;
    }
    for(k = 0; k < k_tc_re[i][j]; ++k){
       KMAXtcre=KMAXtcre+1;
       SynapsesTCRE[KMAXtcre-1][2]=a_tc_re[i][j][k].FF;
    }
  }
for(j = 0; j < Mtc1; ++j)
  for(i = 0; i < Mtc; ++i){
    for(k = 0; k < k_re_tc[i][j]; ++k){
       KMAXretc=KMAXretc+1;
       SynapsesRETC[KMAXretc-1][2]=g_re_tc[i][j][k].FF*tc_cell[0][0].S;;       //SynapsesRETC[KMAXretc-1][1]=g_re_tc[i][j][k].FF;
    }
    for(k = 0; k < k_re_tc[i][j]; ++k){
       KMAXretcB=KMAXretcB+1;
       SynapsesRETCB[KMAXretcB-1][2]=gb[i][j][k].FF;
    }
    for(k = 0; k < k_tc_tc[i][j]; ++k){
       KMAXtctc=KMAXtctc+1;
       SynapsesTCTC[KMAXtctc-1][2]=a_tc_tc[i][j][k].FF;
    }
  }
}

if( (t>(9200-0.5*h)) && (t<(9200+0.5*h)) ){
KMAXrere=0;
KMAXtcre=0;
KMAXretc=0;
KMAXtctc=0;
KMAXretcB=0;
for(j = 0; j < Mre1; ++j)
  for(i = 0; i < Mre; ++i){
    for(k = 0; k < k_re_re[i][j]; ++k){
       KMAXrere=KMAXrere+1;
       SynapsesRERE[KMAXrere-1][3]=g_re_re[i][j][k].FF*re_cell[0][0].S;      //SynapsesRERE[KMAXrere-1][1]=g_re_re[i][j][k].FF;
    }
    for(k = 0; k < k_tc_re[i][j]; ++k){
       KMAXtcre=KMAXtcre+1;
       SynapsesTCRE[KMAXtcre-1][3]=a_tc_re[i][j][k].FF;
    }
  }
for(j = 0; j < Mtc1; ++j)
  for(i = 0; i < Mtc; ++i){
    for(k = 0; k < k_re_tc[i][j]; ++k){
       KMAXretc=KMAXretc+1;
       SynapsesRETC[KMAXretc-1][3]=g_re_tc[i][j][k].FF*tc_cell[0][0].S;;       //SynapsesRETC[KMAXretc-1][1]=g_re_tc[i][j][k].FF;
    }
    for(k = 0; k < k_re_tc[i][j]; ++k){
       KMAXretcB=KMAXretcB+1;
       SynapsesRETCB[KMAXretcB-1][3]=gb[i][j][k].FF;
    }
    for(k = 0; k < k_tc_tc[i][j]; ++k){
       KMAXtctc=KMAXtctc+1;
       SynapsesTCTC[KMAXtctc-1][3]=a_tc_tc[i][j][k].FF;
    }
  }
}

if( (t>(15200-0.5*h)) && (t<(15200+0.5*h)) ){
KMAXrere=0;
KMAXtcre=0;
KMAXretc=0;
KMAXtctc=0;
KMAXretcB=0;
for(j = 0; j < Mre1; ++j)
  for(i = 0; i < Mre; ++i){
    for(k = 0; k < k_re_re[i][j]; ++k){
       KMAXrere=KMAXrere+1;
       SynapsesRERE[KMAXrere-1][4]=g_re_re[i][j][k].FF*re_cell[0][0].S;      //SynapsesRERE[KMAXrere-1][1]=g_re_re[i][j][k].FF;;
    }
    for(k = 0; k < k_tc_re[i][j]; ++k){
       KMAXtcre=KMAXtcre+1;
       SynapsesTCRE[KMAXtcre-1][4]=a_tc_re[i][j][k].FF;
    }
  }
for(j = 0; j < Mtc1; ++j)
  for(i = 0; i < Mtc; ++i){
    for(k = 0; k < k_re_tc[i][j]; ++k){
       KMAXretc=KMAXretc+1;
       SynapsesRETC[KMAXretc-1][4]=g_re_tc[i][j][k].FF*tc_cell[0][0].S;;       //SynapsesRETC[KMAXretc-1][1]=g_re_tc[i][j][k].FF;;
    }
    for(k = 0; k < k_re_tc[i][j]; ++k){
       KMAXretcB=KMAXretcB+1;
       SynapsesRETCB[KMAXretcB-1][4]=gb[i][j][k].FF;
    }
    for(k = 0; k < k_tc_tc[i][j]; ++k){
       KMAXtctc=KMAXtctc+1;
       SynapsesTCTC[KMAXtctc-1][4]=a_tc_tc[i][j][k].FF;
    }
  }
}

//-----------------------------------
     avr=0;
     av=0;
     for(j = 0; j < Mre1; ++j)
       for(i = 0; i < Mre; ++i)
         avr=avr + y_ini[no_re[i][j][0]];
     for(j = 0; j < Mtc1; ++j)
       for(i = 0; i < Mtc; ++i)
         av=av + y_ini[no_tc[i][j][0]];
     av = av/(Mtc1*Mtc);
     avr = avr/(Mre1*Mre);

   if((ii/(ih*1000))*(ih*1000) == ii) {
      printf("\n T= %lf ",t);
//        for(j = 0; j < Mtc1; ++j)
        for(i = 0; i < Mtc; ++i){
        if(I_TC == 1){
          printf("\n TC: ");
          for(k = 0; k < N_TC; ++k) 
             printf("%lf ", y_ini[no_tc[i][0][k]]);}
}
//        for(j = 0; j < Mre1; ++j)
        for(i = 0; i < Mre; ++i){
        if(I_RE == 1){
          printf("\n RE: ");
          for(k = 0; k < N_RE; ++k)   
             printf("%lf ", y_ini[no_re[i][0][k]]);}
}
//        for(j = 0; j < Mc1; ++j)
        for(i = 0; i < Mc; ++i){
        if(I_CX == 1){
          printf("\n CX: ");
          for(k = 0; k < N_CX; ++k) 
             printf("%lf ", cx_cell[i][0].v_SOMA);}
        if(I_IN == 1){
          printf("\n IN: ");
          for(k = 0; k < N_IN; ++k)   
             printf("%lf ", y_ini[no_in[i][j][k]]);}
}
}
   /*  
   if((t > t3D) && ((ii/(ih))*(ih) == ii)){
     for(j = 0; j < Mre1; ++j){
       for(i = 0; i < Mre; ++i){
             if(I_RE == 1) fprintf(f6,"%lf ", y_ini[no_re[i][j][0]]);}
       fprintf(f6,"\n");}

     for(j = 0; j < Mtc1; ++j){
       for(i = 0; i < Mtc; ++i){
             if(I_TC == 1) fprintf(f8,"%lf ", y_ini[no_tc[i][j][0]]); }
       fprintf(f8,"\n");}

     for(j = 0; j < Mc1; ++j){
       for(i = 0; i < Mc; ++i){
//         fprintf(f10,"%lf ", y_ini[no_cx[i][j][0]]); 
         if(I_CX == 1) fprintf(f10,"%lf ", cx_cell[i][j].v_SOMA); 
         if(I_IN == 1) fprintf(f12,"%lf ", in_cell[i][j].v_SOMA); 
         }
       fprintf(f10,"\n");
       fprintf(f12,"\n");
       }
  }
  */
   if((t > ttime) && ((ii/25)*(25) == ii)){
     if(I_RE == 1) fprintf(f7,"%lf %lf %lf ", t, avr,
                             aPRE[0][0].RS*a_ext1.g/re_cell[0][0].S);
			   //                  gb_re_re[0][0][0].G);
     if(I_TC == 1) fprintf(f9,"%lf %lf %lf %lf %lf ", t, av,
                             aP.RS*a_ext1.g/tc_cell[0][0].S,
                             gb[0][0][0].G, g_re_tc[0][0][0].G);
     for(i = 0; i < Mre; ++i){
         if(I_RE == 1) fprintf(f7,"%lf ", y_ini[no_re[i][0][0]]);}
     fprintf(f7,"\n");

     for(i = 0; i < Mtc; ++i){
         if(I_TC == 1) fprintf(f9,"%lf ", y_ini[no_tc[i][0][0]]);}
     fprintf(f9,"\n");
   }

   if((t > ttime) && ((ii/(2))*(2) == ii)){
     if(I_CX == 1) fprintf(f11,"%lf ", t);
     if(I_IN == 1) fprintf(f13,"%lf ", t);
     for(i = 0; i < Mc; ++i){
         if(I_CX == 1) fprintf(f11,"%lf ", cx_cell[i][Mc1/2].v_SOMA);
         if(I_IN == 1) fprintf(f13,"%lf ", in_cell[i][Mc1/2].v_SOMA);
     }
     fprintf(f11,"\n");
     fprintf(f13,"\n");
   }

}

fsRERE = fopen("SYnapceRERE", "w");
fsRETC = fopen("SYnapceRETC", "w");
fsTCRE = fopen("SYnapceTCRE", "w");
fsTCTC = fopen("SYnapceTCTC", "w");
fsRETCB = fopen("SYnapceRETCB", "w"); 

for(k=0; k < KMAXrere; k++){
   for(j=0; j < 5; j++)
      fprintf(fsRERE,"%lf ",SynapsesRERE[k][j]);
   fprintf(fsRERE,"\n ");
} 
for(k=0; k < KMAXretc; k++){
   for(j=0; j < 5; j++)
      fprintf(fsRETC,"%lf ",SynapsesRETC[k][j]);
   fprintf(fsRETC,"\n ");
}
for(k=0; k < KMAXretcB; k++){
   for(j=0; j < 5; j++)
      fprintf(fsRETCB,"%lf ",SynapsesRETCB[k][j]);
   fprintf(fsRETCB,"\n ");
}
for(k=0; k < KMAXtcre; k++){
   for(j=0; j < 5; j++)
      fprintf(fsTCRE,"%lf ",SynapsesTCRE[k][j]);
   fprintf(fsTCRE,"\n ");
}
for(k=0; k < KMAXtctc; k++){
   for(j=0; j < 5; j++)
      fprintf(fsTCTC,"%lf ",SynapsesTCTC[k][j]);
   fprintf(fsTCTC,"\n ");
}

fclose(fsRERE);
fclose(fsRETC);
fclose(fsTCRE);
fclose(fsTCTC);
fclose(fsRETCB);
//--------------------END CALCULATION-------------------------------

//-----------------close ALL files-----------------------------------
  fclose(f2);
  fclose(f22);
  fclose(f6);
  fclose(f7);
  fclose(f8);
  fclose(f9);
  fclose(f10);
  fclose(f11);
  fclose(f12);
  fclose(f13);
  //  fclose(fNoise);
  fclose(fre); 
  printf("\n"); 
}

//+++++++++++ Function to calculate the index of BOUNDARY elements ++++++++++
int b(unsigned btype, unsigned m, int ind) {
      if(btype == 0) {  //------flow boundary conditions:
        if( (ind >= 0) && (ind <= (m-1)) ) return( ind );
        if(ind < 0) return( -ind );
        if(ind > (m-1)) return( 2*m - ind - 2); }

      else if(btype == 1) {  //------periodic boundary conditions: 
        if( (ind >= 0) && (ind <= (m-1)) ) return( ind );
        if(ind < 0) return( ind + m);
        if(ind > (m-1)) return (ind - m); }

      else if(btype == 9) {  //------NO boundary conditions:
        return( ind ); }     
}

//+++++++++++ Function to calculate the right sides for ALL ODE +++++++++++++++
void fun(unsigned neq, double x, double *y_ini, double *f_ini){
double scale;
int i, j, k, i1, j1, ii, jj, nst, kk[Max][Max1];

//========here the MAIN loop to calculate intrinsic conductances===========
//--------(f_ini IS changed, y_ini IS NOT changed)-------------------------
for(i=0; i < Mre; ++i)
for(j=0; j < Mre1; ++j){
  if(I_RE == 1) re_cell[i][j].calc(x, y_ini+no_re[i][j][0], 
                                      f_ini+no_re[i][j][0]); }

for(i=0; i < Mtc; ++i)
for(j=0; j < Mtc1; ++j){
  if(I_TC == 1) tc_cell[i][j].calc(x, y_ini+no_tc[i][j][0], 
                                      f_ini+no_tc[i][j][0]); }

for(i=0; i < Mc; ++i)
for(j=0; j < Mc1; ++j){
  if(I_CX == 1) cx_cell[i][j].calc(x, y_ini+no_cx[i][j][0], 
                                      f_ini+no_cx[i][j][0]);
  if(I_IN == 1) in_cell[i][j].calc(x, y_ini+no_in[i][j][0], 
                                      f_ini+no_in[i][j][0]); }

//========here the MAIN loop to calculate synaptic conductances=============
//--------(f_ini IS changed, y_ini IS NOT changed) -------------------------
for(i = 0; i < Mre; ++i)
for(j = 0; j < Mre1; ++j){

//--------reciprocal GABA-A between RE cells---------------------------------
if(I_RE == 1){
  for(i1 = 0, k = 0; i1 < Mre; ++i1)
  for(j1 = 0; j1 < Mre1; ++j1){

//    scale = sqrt((double) (i1*i1 + j1*j1));
   if(C_RERE[i1][j1][i][j] > 0){
    ii = i1; jj = j1;
//  if( (BOUND != 9) || ( (ii>=0) && (ii<=(Mre-1)) && 
//                                   (jj>=0) && (jj<=(Mre1-1)) ) ){
          if( !( (ii == i) && (jj == j) && (SELFING == 0) ) ){
              g_re_re[i][j][k].calc(C_RERE[i1][j1][i][j]/re_cell[0][0].S, x, y_ini+no_ga_rere[i][j][k][0],    //g_re_re[i][j][k].calc(g_gaba_a*C_RERE[i1][j1][i][j], x, y_ini+no_ga_rere[i][j][k][0],
                      f_ini+no_ga_rere[i][j][k][0],
                      y_ini[no_re[i][j][0]], y_ini[no_re[ii][jj][0]]);
          if(I_GB_RERE == 1){ gb_re_re[i][j][k].calc(g_gaba_b_re_re, x, 
                   y_ini+no_gb_rere[i][j][k][0], f_ini+no_gb_rere[i][j][k][0], 
                   y_ini[no_re[i][j][0]], y_ini[no_re[ii][jj][0]]); }
              ++k; }   // number of synapces of the cell i-j 
// }
  } }
  k_re_re[i][j] = k;
  for(k = 0; k < k_re_re[i][j]; ++k){
    if(k_re_re[i][j] == 0) { }  //NO synapces
    else { 
       g_re_re[i][j][k].I = g_re_re[i][j][k].I / k_re_re[i][j];
       f_ini[no_re[i][j][0]] = f_ini[no_re[i][j][0]] - g_re_re[i][j][k].I; 
      if(I_GB_RERE == 1) {gb_re_re[i][j][k].I = gb_re_re[i][j][k].I / 
                                                k_re_re[i][j];
                          f_ini[no_re[i][j][0]] = f_ini[no_re[i][j][0]] - 
			  gb_re_re[i][j][k].I; }
}
}}

//-----------AMPA from TC to RE cells---------------------------------------
if(I_RE == 1 && I_TC == 1){
  for(i1 = 0, k = 0; i1 < Mtc; ++i1)
  for(j1 = 0; j1 < Mtc1; ++j1){
//    scale = sqrt((double) (i1*i1 + j1*j1));
//    if(scale <= MS_TC_RE_MAX){
   if(C_TCRE[i1][j1][i][j] > 0){
    ii = i1; jj = j1;
//  if( (BOUND != 9) || ( (ii>=0) && (ii<=(Mtc-1)) && 
//                                   (jj>=0) && (jj<=(Mtc1-1)) ) ){
        a_tc_re[i][j][k].calc(g_ampa, x, y_ini[no_re[i][j][0]], 
                           y_ini[no_tc[ii][jj][0]]);
        ++k;  // number of synapces of the cell i-j 
// }
  } }
  k_tc_re[i][j] = k;
  for(k = 0; k < k_tc_re[i][j]; ++k){
    if(k_tc_re[i][j] == 0) { }  //NO synapces
    else {
       a_tc_re[i][j][k].I = a_tc_re[i][j][k].I / k_tc_re[i][j];
       f_ini[no_re[i][j][0]] = f_ini[no_re[i][j][0]] - a_tc_re[i][j][k].I; }
}}

//-------------AMPA from CX to RE cells-------------------------------------
if(I_CX == 1 && I_RE == 1){
  for(i1 = -MS_CX_TC, k = 0; i1 <= MS_CX_TC; ++i1)
  for(j1 = -MS_CX_TC1; j1 <= MS_CX_TC1; ++j1){
    scale = sqrt((double) (i1*i1 + j1*j1));
    if(scale <= MS_CX_TC_MAX){
      ii = i + i1; jj = j + j1;
      if( (BOUND != 9) || ( (ii>=0) && (ii<=(Mc-1)) && 
                                       (jj>=0) && (jj<=(Mc1-1)) ) ){
        a_cx_re[i][j][k].calc(g_ampa_cx_re, x, y_ini[no_re[i][j][0]], 
                           cx_cell[b(BOUND,Mc,ii)][b(BOUND,Mc1,jj)].v_SOMA);
        ++k; } // number of synapces of the cell i-j 
  } }
  kk[i][j] = k;
  for(k = 0; k < kk[i][j]; ++k){
    if(kk[i][j] == 0) { }  //NO synapces
    else {
       a_cx_re[i][j][k].I = a_cx_re[i][j][k].I / kk[i][j];
       f_ini[no_re[i][j][0]] = f_ini[no_re[i][j][0]] - a_cx_re[i][j][k].I; }
}}
}


for(i = 0; i < Mtc; ++i)
for(j = 0; j < Mtc1; ++j){

//--------reciprocal AMPA between TC cells----------------------------------
 if(I_TC == 1){
  for(i1 = 0, k = 0; i1 < Mtc; ++i1)
  for(j1 = 0; j1 < Mtc1; ++j1){

//    scale = sqrt((double) (i1*i1 + j1*j1));
//    if(scale <= MS_TC_TC_MAX){
   if(C_TCTC[i1][j1][i][j] > 0){
    ii = i1; jj = j1;
//  if( (BOUND != 9) || ( (ii>=0) && (ii<=(Mtc-1)) && 
//                                     (jj>=0) && (jj<=(Mtc1-1)) ) ){
          if( !( (ii == i) && (jj == j) && (SELFING == 0) ) ){
              a_tc_tc[i][j][k].calc(g_ampa1, x,  y_ini[no_tc[i][j][0]], 
                               y_ini[no_tc[ii][jj][0]]);
              ++k; }   // number of synapces of the cell i-j 
// }
  } }
  k_tc_tc[i][j] = k;
  for(k = 0; k < k_tc_tc[i][j]; ++k){
    if(k_tc_tc[i][j] == 0) { }  //NO synapces
    else { 
       a_tc_tc[i][j][k].I = a_tc_tc[i][j][k].I / k_tc_tc[i][j];
       f_ini[no_tc[i][j][0]] = f_ini[no_tc[i][j][0]] - a_tc_tc[i][j][k].I; }

}}

//---------GABA-A and GABA-B from RE to TC cells------------------------------
if(I_RE == 1 && I_TC == 1){
  for(i1 = 0, k = 0; i1 < Mre; ++i1)
  for(j1 = 0; j1 < Mre1; ++j1){
//    scale = sqrt((double) (i1*i1 + j1*j1));
//    if(scale <= MS_RE_TC_MAX){
   if(C_RETC[i1][j1][i][j] > 0){
    ii = i1; jj = j1;
//  if( (BOUND != 9) || ( (ii>=0) && (ii<=(Mre-1)) && 
//                                     (jj>=0) && (jj<=(Mre1-1)) ) ){
       g_re_tc[i][j][k].calc(C_RETC[i1][j1][i][j]/tc_cell[0][0].S, x, y_ini+no_ga_retc[i][j][k][0],    //g_re_tc[i][j][k].calc(g_gaba_a1*C_RETC[i1][j1][i][j], x, y_ini+no_ga_retc[i][j][k][0],
                      f_ini+no_ga_retc[i][j][k][0],
                      y_ini[no_tc[i][j][0]], y_ini[no_re[ii][jj][0]]);
       if(I_GB == 1){ gb[i][j][k].calc(g_gaba_b, x, y_ini+no_g[i][j][k][0],
                           f_ini+no_g[i][j][k][0], y_ini[no_tc[i][j][0]], 
                             y_ini[no_re[ii][jj][0]]); }
         ++k;    // number of synapces of the cell i-j
// }
  } }
  k_re_tc[i][j] = k;
  for(k = 0; k < k_re_tc[i][j]; ++k){
    if(k_re_tc[i][j] == 0) { }  //NO synapces
    else {
       g_re_tc[i][j][k].I = g_re_tc[i][j][k].I / k_re_tc[i][j];
       if(I_GB == 1) gb[i][j][k].I = gb[i][j][k].I / k_re_tc[i][j];
       f_ini[no_tc[i][j][0]] = f_ini[no_tc[i][j][0]] - g_re_tc[i][j][k].I;
       if(I_GB == 1) f_ini[no_tc[i][j][0]] = f_ini[no_tc[i][j][0]] - 
                                                          gb[i][j][k].I; 
}
}}

//------------AMPA from CX to TC cells-------------------------------------
if(I_CX == 1 && I_TC == 1){
  for(i1 = -MS_CX_TC, k = 0; i1 <= MS_CX_TC; ++i1)
  for(j1 = -MS_CX_TC1; j1 <= MS_CX_TC1; ++j1){
    scale = sqrt((double) (i1*i1 + j1*j1));
    if(scale <= MS_CX_TC_MAX){
      ii = i + i1; jj = j + j1;
      if( (BOUND != 9) || ( (ii>=0) && (ii<=(Mc-1)) && 
                                       (jj>=0) && (jj<=(Mc1-1)) ) ){
        a_cx_tc[i][j][k].calc(g_ampa_cx_tc, x, y_ini[no_tc[i][j][0]], 
                           cx_cell[b(BOUND,Mc,ii)][b(BOUND,Mc1,jj)].v_SOMA);
        ++k; } // number of synapces of the cell i-j 
  } }
  kk[i][j] = k;
  for(k = 0; k < kk[i][j]; ++k){
    if(kk[i][j] == 0) { }  //NO synapces
    else {
       a_cx_tc[i][j][k].I = a_cx_tc[i][j][k].I / kk[i][j];
       f_ini[no_tc[i][j][0]] = f_ini[no_tc[i][j][0]] - a_cx_tc[i][j][k].I; }
}}
}

for(i = 0; i < Mc; ++i)
for(j = 0; j < Mc1; ++j){

//-----------reciprocal AMPA 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){
    scale = sqrt((double) (i1*i1 + j1*j1));
    if(scale <= MS_CX_CX_MAX){
      ii = i + i1; jj = j + j1;
      if( (BOUND != 9) || ( (ii>=0) && (ii<=(Mc-1)) && 
                                       (jj>=0) && (jj<=(Mc1-1)) ) ){
        if( !( ((i1*i1 + j1*j1) == 0) && (SELFINGcx == 0) ) ){
        a_cx_cx[i][j][k].calc(g_ampa_cx_cx, x, y_ini[no_cx[i][j][0]], 
                           cx_cell[b(BOUND,Mc,ii)][b(BOUND,Mc1,jj)].v_SOMA);
        ++k; }} // number of synapces of the cell i-j 
  } }
  kk[i][j] = k;
  for(k = 0; k < kk[i][j]; ++k){
    if(kk[i][j] == 0) { }  //NO synapces
    else {
       a_cx_cx[i][j][k].I = a_cx_cx[i][j][k].I / kk[i][j];
       f_ini[no_cx[i][j][0]] = f_ini[no_cx[i][j][0]] - a_cx_cx[i][j][k].I; }
}}

//-------------AMPA 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){
    scale = sqrt((double) (i1*i1 + j1*j1));
    if(scale <= MS_CX_IN_MAX){
      ii = i + i1; jj = j + j1;
      if( (BOUND != 9) || ( (ii>=0) && (ii<=(Mc-1)) && 
                                       (jj>=0) && (jj<=(Mc1-1)) ) ){
        a_cx_in[i][j][k].calc(g_ampa_cx_in, x, y_ini[no_in[i][j][0]], 
                           cx_cell[b(BOUND,Mc,ii)][b(BOUND,Mc1,jj)].v_SOMA);
        ++k; } // number of synapces of the cell i-j 
  } }
  kk[i][j] = k;
  for(k = 0; k < kk[i][j]; ++k){
    if(kk[i][j] == 0) { }  //NO synapces
    else {
       a_cx_in[i][j][k].I = a_cx_in[i][j][k].I / kk[i][j];
       f_ini[no_in[i][j][0]] = f_ini[no_in[i][j][0]] - a_cx_in[i][j][k].I; }
}}

//--------------GABA-A 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){
    scale = sqrt((double) (i1*i1 + j1*j1));
    if(scale <= MS_IN_CX_MAX){
      ii = i + i1; jj = j + j1;
      if( (BOUND != 9) || ( (ii>=0) && (ii<=(Mc-1)) && 
                                       (jj>=0) && (jj<=(Mc1-1)) ) ){
        ga_in_cx[i][j][k].calc(g_gaba_a_in_cx, x, y_ini[no_cx[i][j][0]], 
                           in_cell[b(BOUND,Mc,ii)][b(BOUND,Mc1,jj)].v_SOMA);
//                           y_ini[no_in[b(BOUND,Mc,ii)][b(BOUND,Mc1,jj)][0]]);
        ++k; } // number of synapces of the cell i-j 
  } }
  kk[i][j] = k;
  for(k = 0; k < kk[i][j]; ++k){
    if(kk[i][j] == 0) { }  //NO synapces
    else {
       ga_in_cx[i][j][k].I = ga_in_cx[i][j][k].I / kk[i][j];
       f_ini[no_cx[i][j][0]] = f_ini[no_cx[i][j][0]] - ga_in_cx[i][j][k].I; }
}}

//--------------GABA-B 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){
    scale = sqrt((double) (i1*i1 + j1*j1));
    if(scale <= MS_IN_CX_MAX){
      ii = i + i1; jj = j + j1;
      if( (BOUND != 9) || ( (ii>=0) && (ii<=(Mc-1)) && 
                                       (jj>=0) && (jj<=(Mc1-1)) ) ){
        gb_in_cx[i][j][k].calc(g_gaba_b_in_cx, x, y_ini+no_gcx[i][j][k][0],
                           f_ini+no_gcx[i][j][k][0], y_ini[no_cx[i][j][0]], 
                           in_cell[b(BOUND,Mc,ii)][b(BOUND,Mc1,jj)].v_SOMA);
//                           y_ini[no_in[b(BOUND,Mc,ii)][b(BOUND,Mc1,jj)][0]]);
        ++k; } // number of synapces of the cell i-j 
  } }
  kk[i][j] = k;
  for(k = 0; k < kk[i][j]; ++k){
    if(kk[i][j] == 0) { }  //NO synapces
    else {
       gb_in_cx[i][j][k].I = gb_in_cx[i][j][k].I / kk[i][j];
       f_ini[no_cx[i][j][0]] = f_ini[no_cx[i][j][0]] - gb_in_cx[i][j][k].I; }
}}

//--------------AMPA from TC to CX cells-------------------------------------
if(I_CX == 1 && I_TC == 1){
  for(i1 = -MS_TC_CX, k = 0; i1 <= MS_TC_CX; ++i1)
  for(j1 = -MS_TC_CX1; j1 <= MS_TC_CX1; ++j1){
    scale = sqrt((double) (i1*i1 + j1*j1));
    if(scale <= MS_TC_CX_MAX){
      ii = i + i1; jj = j + j1;
      if( (BOUND != 9) || ( (ii>=0) && (ii<=(Mtc-1)) && 
                                       (jj>=0) && (jj<=(Mtc1-1)) ) ){        
        a_tc_cx[i][j][k].calc(g_ampa_tc_cx, x, y_ini[no_cx[i][j][0]], 
                           y_ini[no_tc[b(BOUND,Mtc,ii)][b(BOUND,Mtc1,jj)][0]]);
        ++k; } // number of synapces of the cell i-j 
  } }
  kk[i][j] = k;
  for(k = 0; k < kk[i][j]; ++k){
    if(kk[i][j] == 0) { }  //NO synapces
    else {
       a_tc_cx[i][j][k].I = a_tc_cx[i][j][k].I / kk[i][j];
       f_ini[no_cx[i][j][0]] = f_ini[no_cx[i][j][0]] - a_tc_cx[i][j][k].I; }
}}
//-------------AMPA from TC to IN cells--------------------------------------
if(I_IN == 1 && I_TC == 1){
  for(i1 = -MS_TC_IN, k = 0; i1 <= MS_TC_IN; ++i1)
  for(j1 = -MS_TC_IN1; j1 <= MS_TC_IN1; ++j1){
    scale = sqrt((double) (i1*i1 + j1*j1));
    if(scale <= MS_TC_IN_MAX+0.0001){
      ii = i + i1; jj = j + j1;
      if( (BOUND != 9) || ( (ii>=0) && (ii<=(Mtc-1)) && 
                                       (jj>=0) && (jj<=(Mtc1-1)) ) ){
       a_tc_in[i][j][k].calc(g_ampa_tc_in, x, y_ini[no_in[i][j][0]], 
                           y_ini[no_tc[b(BOUND,Mtc,ii)][b(BOUND,Mtc1,jj)][0]]);
        ++k; } // number of synapces of the cell i-j 
  } }
  kk[i][j] = k;
  for(k = 0; k < kk[i][j]; ++k){
    if(kk[i][j] == 0) { }  //NO synapces
    else {
       a_tc_in[i][j][k].I = a_tc_in[i][j][k].I / kk[i][j];
       f_ini[no_in[i][j][0]] = f_ini[no_in[i][j][0]] - a_tc_in[i][j][k].I; }
}}

}
//=============ENF of MAIN loop==================================================

//------------------expernal stimulation---------------------------------
nst = 0;   //type of external stimulation

//-----------train of shocks----------------------------------------
if(nst == 0){ 

  for(i = 0; i < Mtc; ++i)
    for(j = 0; j < Mtc1; ++j){
      a_ext3.calc(g_ext3, x);
      if(I_TC == 1) f_ini[no_tc[i][j][0]] = f_ini[no_tc[i][j][0]] - 
	                     y_ini[no_tc[i][j][0]] * a_ext3.g/tc_cell[0][0].S;
    }

  for(i = 0; i < Mre; ++i)
    for(j = 0; j < Mre1; ++j){
      a_ext4.calc(g_ext4, x);
      if(I_RE == 1) f_ini[no_re[i][j][0]] = f_ini[no_re[i][j][0]] - 
	                     y_ini[no_re[i][j][0]] * a_ext4.g/re_cell[0][0].S;
    }

//      if(g_ext1 > 1.e-10){
//      fscanf(fNoise,"%lf %lf",&tt, &Noise);
  a_ext1.Noise = Noise;
  a_ext2.Noise = Noise;
      a_ext1.calc(g_ext1, x);
      a_ext2.calc(g_ext2, x);

    aP.calc(x);


      
//    if(I_TC == 1) a_ext1.calc(g_ext_tc, x);
//    if(I_RE == 1) a_ext2.calc(g_ext_re, x);
//    if(I_CX == 1) a_ext3.calc(g_ext_cx, x);
//    if(I_IN == 1) a_ext4.calc(g_ext_in, x);

      for(i = 0; i < Mc; ++i)
      for(j = 0; j < Mc1; ++j){
         if(I_IN == 1) f_ini[no_in[i][j][0]] = f_ini[no_in[i][j][0]] - 
                       y_ini[no_in[i][j][0]] * a_ext4.g; 
         if(I_CX == 1) f_ini[no_cx[i][j][0]] = f_ini[no_cx[i][j][0]] - 
                       y_ini[no_cx[i][j][0]] * a_ext3.g; 
      }
      for(i = 0; i < Mtc; ++i)
      for(j = 0; j < Mtc1; ++j){
        if(C_TC[i][j] > 0)
           if(I_TC == 1) f_ini[no_tc[i][j][0]] = f_ini[no_tc[i][j][0]] - 
	             C_TC[i][j]*y_ini[no_tc[i][j][0]] * aP.RS * a_ext1.g/tc_cell[0][0].S;
           }

      for(i = 0; i < Mre; ++i)
      for(j = 0; j < Mre1; ++j){
         aPRE[i][j].calc(x);
         if(C_RE[i][j] > 0){
           if(I_RE == 1) f_ini[no_re[i][j][0]] = f_ini[no_re[i][j][0]] - 
	     C_RE[i][j]*y_ini[no_re[i][j][0]] * aPRE[i][j].RS * a_ext2.g/re_cell[0][0].S;}
	   } 
      //     }  
} 
















//-----------ONE shock--------------------------------------------- 
else if(nst == 1) {
  if( ((x > 2000) && (x < 2050)) || ((x > 9000) && (x < 9050)) ) {
     a_ext1.calc(g_ext_tc, x);
  
        if(I_TC == 1) f_ini[no_tc[0][0][0]] = f_ini[no_tc[0][0][0]] -
           a_ext1.g * y_ini[no_tc[0][0][0]];
  }
}

//-----------temporal stuff------------------------------------------
//if(x > 1000) cx_cell[0][0].I_Stim1 = 0.375; //0.5 would correspond to 100 for soma
//if(x > 2000) cx_cell[0][0].I_Stim1 = 0.;

//if(x > 1000) cx_cell[0][0].I_Stim2 = 70;
//if(x > 2000) cx_cell[0][0].I_Stim2 = 0;

//if(x > 1000) cx_cell[0][0].I_Stim2 = 600;
//if(x > 1020) cx_cell[0][0].I_Stim2 = 0;

//if((x > 1000) && (x < 2000)) f_ini[no_cx[0][0][0]] = f_ini[no_cx[0][0][0]] + 
//                                                                      (0.666666);
//if((x > 3000) && (x < 4000)) cx_cell[0][0].v_SOMA = cx_cell[0][0].v_SOMA + (1.);
//if((x > 1900) && (x < 2000)) re_cell[0][0].f[0] = re_cell[0][0].f[0] - (1.3);
//if((x > 1000) && (x < 1200)) y_ini[no_in[0][0][0]] = y_ini[no_in[0][0][0]] + (0.5);
//if((x > 7000) && (x < 7200)) re_cell[0][0].f[0] = re_cell[0][0].f[0] - (0.5);
//if((x > 0) && (x < 100)) tc_cell[0][0].f[0] = tc_cell[0][0].f[0] - (0.5);
//if((x > 2000) && (x < 2100)) tc_cell[0][0].f[0] = tc_cell[0][0].f[0] + (0.3);
//if((x > 1200) && (x < 1250)) tc_cell[0][0].f[0] = tc_cell[0][0].f[0] - (0.3);
//if((x > 1440) && (x < 1500)) tc_cell[0][0].f[0] = tc_cell[0][0].f[0] - (0.3);
//if((x > 1700) && (x < 1800)) tc_cell[0][0].f[0] = tc_cell[0][0].f[0] - (0.3);
//if((x > 10000) && (x < 11000)) tc_cell[0][0].f[0] = tc_cell[0][0].f[0] - (0.2);
//if((x > 15000) && (x < 15500)) tc_cell[0][0].f[0] = tc_cell[0][0].f[0] + (1.7);
//if((x > 2500) && (x < 3000)) tc_cell[0][0].f[0] = tc_cell[0][0].f[0] + (0.2);
}

/*
void solout (long nr, double xold, double x, double* y, unsigned n, int* irtrn) 
{ 
  static double xout, xout1; 
  char format99[] = "x=%f  y=%12.10f %12.10f  nstep=%li\r\n"; 
  int i, j, k;

  if (nr == 1) {  
//    printf ( "x=%f  y=%12.10f %12.10f  nstep=%li\r\n", x, y[0], y[1], nr-1); 

    fprintf(f6,"%lf ", x);
    fprintf(f8,"%lf ", x);
    for(j = 0; j < M1; ++j){
       for(i = 0; i < M; ++i){
         fprintf(f6,"%lf ", y[no_re[i][j]]);
         fprintf(f8,"%lf ", y[no_tc[i][j]]); }
    fprintf(f6,"\n");
    fprintf(f8,"\n");}
       
    xout = x + 0.2; 
  } 
  else {
    while (x >= xout1) {
      printf("\n t=%lf", xout1);
      xout1 += 100.0; } 

    while (x >= xout) { 
//      printf (format99, xout, contd8(0,xout), contd8(1,xout), nr-1); 

      fprintf(f6,"%lf ", xout);
      fprintf(f8,"%lf ", xout);
      for(j = 0; j < M1; ++j){
         for(i = 0; i < M; ++i){
           fprintf(f6,"%lf ", contd8(no_re[i][j],xout) );
           fprintf(f8,"%lf ", contd8(no_tc[i][j],xout) ); }
      fprintf(f6,"\n");
      fprintf(f8,"\n");}

      xout += 0.2; 
    } 
  }
} */














Loading data, please wait...