main.cpp

#include <iostream.h>       // for output
#include<fstream.h>         // for file output
#include "sstream.h"        // for string output
#include <string>


#include "CmprtmntRk4.h"    // for cylindrical compartment		    
#include "CurrentRk4.h"

#include "Chanlist.h"       // my library of ion channels  
#include "param.h"          // file for all parameters 

#include "Injector.h"		// for current injection
#include "Input.h"   
#include "Noise.h"          // Cell noise generator
#include "NoiseSource.h"    // Header file for all kind of noise variable

#include <math.h> 
#include<stdlib.h>
#include<time.h>          

#ifdef macintosh
#include <console.h>		// for ccommand() function
#endif



// Reading files
const int read = 0;

//Cells
Compartment cell[NetSize][NetSize];
Input *input[NetSize][NetSize];       // Their inputs

// Connectivity

AlphaSynS *SynS[NetSize][NetSize][NetSize][NetSize];     
real Ji[NetSize][NetSize][NetSize][NetSize];              //synaptic strength (Inhibitory)
    

AlphaSyn *SynE[NetSize][NetSize][NetSize][NetSize];        
real Je[NetSize][NetSize][NetSize][NetSize];              //synaptic strength (Exitatory)

 
// Output files

const int M1=0;                                           //Index of the cells that are recorded
const int K1=0; 
const int M2=0;
const int K2=1;

real LFP1[Nstep];                                         //LFP traces
real LFP2[Nstep];
real LFP3[Nstep];
real LFP4[Nstep];
real LFP5[Nstep];
int counterLFP3[NetSize][NetSize];

real Raster[NetSize*NetSize][1000];                        // raster of spikes
int gRaster[Nstep];                                       // global raster
int SynRaster[Nstep];
int Rcount[Ncell]; 
real tableV[3][Nstep];                                    //For ouputs  
real Vmem[NetSize][NetSize];                             
int count = 0;                                            //Counter for outputs

const char DELIM = '\t';		                          //field delimiter for output

char a[30]="réseau";                                      //Generic name of the output files
char b[100]="réseau.txt";                                 //Char variable for output files

real variable[11];
#include "LFPfiltDecl.h"

/********************* Method for creating mitral cells *********************/
 

void MakeMitral(Compartment *mitral){
    
             
	mitral->SetV( EREST );		// start at rest
    mitral->Gm=Area/RM;
    mitral->Cm=Area*CM;
    mitral->EGm=ELEAK*Area/RM;
    mitral->MaxGI=gSingI*Area;
    mitral->P0=P0;
    mitral->tau1=taurise;
    mitral->tau2=taudecay;
    mitral->EI=Ein;
    
    ChanStd *NaChan = new ChanStd( mitral, gNa*(1+Var()*DeltaG)*Area);
    //SetUpNaBhalaChannelTres( *NaChan,Slop );
	SetUpNaBhalaChannel( *NaChan );
    
    ChanTM  *KfstChan = new ChanTM( mitral, gK*(1+Var()*DeltaG)*Area );
	SetUpKfastChannel( *KfstChan );
    
    ChanTM *NaPChan = new ChanTM( mitral, gNaP*(1+Var()*DeltaG)*Area );
    SetUpNaPWangXJChannel( *NaPChan );
     
    ChanTM *Ks1Chan = new ChanTM( mitral, gKs*(1+Var()*DeltaG)*Area );
    SetUpKs1WangXJChannel( *Ks1Chan );
     
    ChanTM *AChan = new ChanTM( mitral, gA*(1+Var()*DeltaG)*Area );
    SetUpAWangChannel( *AChan );
   
    Noise *noise = new Noise( mitral );
    noise->SetAmplitude(0.008*Area); 
    
}

void SetMitral(Compartment *mitral){
    
             
	mitral->SetV( EREST );		// start at rest
    mitral->Gm=Area/RM;
    mitral->Cm=Area*CM;
    mitral->EGm=ELEAK*Area/RM;
    mitral->MaxGI=gSingI*Area;
    mitral->P0=P0;
    mitral->tau1=taurise;
    mitral->tau2=taudecay;
    mitral->EI=Ein;
 }

/**************************************************************************/
/*************************     Main Program     ***************************/
/**************************************************************************/


main( int argc, char **argv )
{
/************************* Initiates random numbers *********************/
    
    srand(time(NULL));

/*** give Mac Symantec C++ users a chance to redirect input & output  ***/
	#ifdef macintosh
		ccommand( &argv );
	#endif

/************* Create the cells; add their channels and their current injection ***************/   
for(int m1=0; m1<NetSize; m1+=1){      
     for(int k1=0; k1<NetSize; k1+=1){  
         MakeMitral(&cell[m1][k1]); 
     }
}    


/********************************  Preparing Results output file  ***********************/
 	stringstream ss;
	ss<<"Results"<<".txt";
	ss>>b;
    ofstream Results(b);    
    
    Results<<"fLFP"<<DELIM<<"LFPOscIndex"<<DELIM;
    Results<<"index1"<<DELIM<<"1/sqrt(N)"<<DELIM<<"phase"<<DELIM<<"N"<<DELIM;     
    Results<<"mPsth"<<DELIM<<"index2"<<DELIM;
    Results<<"vLFP1"<<DELIM<<"vLFP2"<<DELIM;
    Results<<"F0"<<DELIM<<"index3"<<DELIM;
    Results<<"tau"<<DELIM<<"Minn1"<<DELIM;
    Results<<"F"<<DELIM<<"tau"<<DELIM<<"Minn2"<<DELIM;         
    Results<<"F"<<DELIM<<"tau"<<DELIM<<"Minn3"<<endl;                     
                                                	
/********************************  Preparing Simulations and running...  ***********************/

variable[0]=0;
variable[1]=.1;
variable[2]=.5;
variable[3]=1;
variable[4]=2;
variable[5]=4;
variable[6]=6;
variable[7]=8;



for(int f=0;f<8; f++){ 
for(int o=0; o<10;o++){

gI=0.03;
gSingI=0.007*variable[f];
SgI=gI*10;
gE=0;

for(int m1=0; m1<NetSize; m1+=1){      
     for(int k1=0; k1<NetSize; k1+=1){  
         SetMitral(&cell[m1][k1]); 
     }
}  

/*********************** defining simulation name  *************************/
stringstream s0;
   s0<<"Res"<<"I"<<100*gI<<"SI"<<100*gSingI;
   s0>>a;
   
/***********************  Builds the Network   ************************/
#include "Construction.h" 
/***********************  A file for all parameters  ************************/  
#include "PropOut.h"

	
   
   
   
   /*********************************************************************/	
   /*****             Runs the simulation for all cells             *****/
   /*****                                                           *****/
   /*********************************************************************/


    
             
// Initializes : this is question of life and death. Many of the steppers  
// do not work without initialization
          for(int m1=0; m1<NetSize; m1+=1){
              for(int k1=0; k1<NetSize; k1++){        
                  counterLFP3[m1][k1]=500;
              }
          }
      
          gStepmaster.doStepinit(DT);  
          count=0;
          
          for(int k1=0; k1<TMAX/DT; k1++){
               LFP1[k1]=0;
               LFP2[k1]=0;
               LFP3[k1]=0;
               gRaster[k1]=0;
               SynRaster[k1]=0;
   	      }
   	      
   	      for(int k1=0; k1<NetSize*NetSize;k1++){
               Rcount[k1]=0;
               for(int m=0; m<100; m++){
                   Raster[k1][m]=0;              
               }  	      
   	      }
   	      
//Computes each time step DT record every DT. 
//This is the core of the simulation !!!

	for (real T=0; T<TMAX; T += DT) {
	
		  gStepmaster.doStepk1(DT);
          gStepmaster.doStepk2(DT);
          gStepmaster.doStepk3(DT);
          gStepmaster.doStepk4(DT);
          gStepmaster.StepAll(DT);          
      
        tableV[0][count]= 1000*cell[M1][K1].GetV();
        tableV[1][count]= 1000*cell[M2][K2].GetV();
 
        
  //Analysis in the step
             
        
        for(int m1=0; m1<NetSize; m1++){           // Summate all potentials for the LFP
           for(int k1=0; k1<NetSize; k1++){
             LFP1[count]+=10*(cell[m1][k1].GetV()+0.07);
             counterLFP3[m1][k1]++;
                      
             if(cell[m1][k1].GetV()>-0.054 && Vmem[m1][k1]<-0.054){
               counterLFP3[m1][k1]=0;   
             }
             
             if(counterLFP3[m1][k1]<500){      // for 20ms
               LFP3[count]+=10*0.014;  
             }
             else{       
               LFP3[count]+=10*(cell[m1][k1].GetV()+0.07);   
             } 
                
             
             if(cell[m1][k1].Memory[0]>0){
                Raster[k1+NetSize*m1][Rcount[k1+NetSize*m1]]=count*DT;
                gRaster[count]++;
                Rcount[k1+NetSize*m1]++;
             }
             
             Vmem[m1][k1]=cell[m1][k1].GetV();
                   
                                        
           }
        }
        count+=1;
    }
    
        
/******************   Writes all recorded traces in a  file    ****************/
    stringstream s4;	
    s4<<a<<".txt";	
	s4>>b;
	ofstream exfile(b);
	
#include "LFPfilter.h"
	
	for(int j=0; j<TMAX/DT; j+=1){
	    exfile << tableV[0][j] << DELIM
        <<tableV[1][j]<<DELIM
        <<LFP2[j]<<DELIM
        <<LFP1[j]<<DELIM
        <<gRaster[j]<<DELIM
        <<LFPpic1[j]//LFP3[j]
        <<endl;
	}
	
	exfile.close();
	
	
	stringstream s5;
	s5<<"Rast"<<a<<".txt";
	s5>>b;
	
	ofstream exfileR(b);
	
	
	for(int j=0; j<Ncell; j+=1){
	    for(int k=0; k<Rcount[j]; k++){
	       exfileR<<Raster[j][k]<<DELIM;
	    }
        exfileR<<endl;
	}
	
	exfileR.close();

    count=0;
    for(int k1=0; k1<Ncell;k1++){
               Rcount[k1]=0;
               for(int m=0; m<100; m++){
                   Raster[k1][m]=0;              
               }  	      
   	      }
   	      
   
   
/************************** Clear objects *********************************/ 

#include "DestructClose.h"

	
}//o
}//f


Results.close();

return 1;
}