Short Term Depression, Presynaptic Inhib., Neuron Diversity Roles in Antennal Lobe (Wei & Lo 2020)

 Download zip file 
Help downloading and running models
Accession:238959

Reference:
1 . Wei KK, Lo CC (2020) Short Term Depression, Presynaptic Inhibition and Local Neuron Diversity Play Key Functional Roles in the Insect Antennal Lobe Journal of Computational Neuroscience, accepted
Model Information (Click on a link to find other models with that property)
Model Type:
Brain Region(s)/Organism:
Cell Type(s): Abstract integrate-and-fire leaky neuron;
Channel(s):
Gap Junctions:
Receptor(s): Cholinergic Receptors; GabaA; GabaB; Glutamate; NMDA;
Gene(s):
Transmitter(s): Acetylcholine; Gaba; Glutamate;
Simulation Environment: C or C++ program;
Model Concept(s): Short-term Synaptic Plasticity;
Implementer(s):
Search NeuronDB for information about:  GabaA; GabaB; NMDA; Glutamate; Cholinergic Receptors; Acetylcholine; Gaba; Glutamate;
/
WeiLo2020
readme.txt
cl_realsimu7e12.c
network.conf
network.pro
rando2.h
                            
// cl_realsimu7e12.c:
// Based on cl_realsimu7e10.c. Here I modify the presynaptic inhibition mechanism so that it affects the vesicle release
//
// Notes:
// [Oct 15, 14] ver 7e12: modifying the mechanism of presynaptic inhibition. Now the presynaptic inhibition also affects
//              short-term plasticity by reducing the amount of vesicle release.
// [Oct 12, 13] ver 7e11: adding new parameters "RefractT" for neuron. The parameter indicates the refractory period
//              for each neuron in the unit of time step.
// [Aug 26, 10] ver 7e10: [Ca2+]^n, n=3.5
// [Jun 22, 10] ver 7e9: add pre-synaptic inhibition
// [May 7, 10] ver 7e8: remove voltage dependence in NMDA synaptic efficacy. Change recpetor naming:
//              NMDA->GABAB, GABA->GABAA, AMPA->Ach
// [Mar 13, 10] ver 7e7: add FreqExtDecay as a changable event to allow decaying external frequency
// [Mar 11, 10] ver 7e6: increase maximum number of target population and decrease the
//                      axonal delay for spikes to 1ms
// [Mar 11, 10] ver 7e6: change the connectivity to one-to-one
// [Apr 15, 08] make output file names changable
// [Apr 5, 05]  change the algorithm for the noise generator.
// [Mar 3, 05]  Add FreqExtSD into the structure PopDescr to describe
//              the standard deviation of the time variation of FreqExt.
//              Note that due to the design of the program structure,
//              FreqExtSD is incompatible with ExtEffSD. If FreqExtSD!=0
//              ExtEffSD will not function. All external efficacies will
//              be set to MeanExtEff.
// [Mar 2, 05]  change the outward rectifying K channel to be time-
//              and voltage-dependent, just as described in Compte-A-2003a
// [Feb 28, 05] make the synaptic efficacy changeable throug .pro file:
//              modified: structure EventDescr, function ParseProtocol()
//              and function HandleEvent().
// [Feb 03, 05] add distribution width to the synaptic efficacy
// [Oct 07, 04] Change max event number from 100 to 200 (line 255)
// [Oct 06, 04] Change the short-term facilitation and short-term
//              depression from cell-wide properties to synapses properties.
// [Oct 06, 04] Change synaptic delay from 5ms to 2ms
// [Oct 05, 04] Add short-term facilitation
// [Sep 22, 04] change NumberofTraces from 2 to 5 (line 1197)
// [Sep 22, 04] Add potassium outward rectifying conductances.
// [Sep 21, 04] fix a bug in DescribeNetwork() relating to the initialization
//              of the population parameters. This bug makes the program initialize
//              the short-term-depression related parameters only for the first
//              population and neglect the rest populations.
// [Sep 21, 04] add potassium inward rectifying conductances (ADR).
//
// [Jun 08, 04] modyfying SaveSpikes(). Use the sliding window instead of the
//              exponential filter.
// [Jun 07, 04] adding code to initialize parameters for short-term depression
//              in case that they are not givien in the configuration file.
// [Jun 02, 04] adding short-term depression
// [Mar 04, 04] adding spike-frequency adaptation from the paper Wang-X-1999a
//
// by Chung-Chuan Lo
//
//###########################################
// Original notes from realsimu8.c
//
// realsimu8 - Ver 0.8
//
// 23-9-03 Number of synapses in GenerateNetwork: tested
// 28-9-03 Parser for the description introduced
// 29-9-03 Protocol parser added
// 30-9-03 Mg block added. Vaux added
// 1-10-03 NMDA saturation
// 2-10-03 bug (use of Tableofspikes) fixed. Sourceneuron introduced
// 8-10-03 bug fixed: two TimeLastemitted spikes are necessary, otherwise the decay for NMDA s variable is always the delay
// 13-10-03 bug fixed: DeltaLS should be \sum_k \delta_k (1-s_k)alpha w_k and not the sum of all s_k
// 15-10-03 successful test for XJ model
// 15-10-03 multitrial version (ver 0.8)

/* How to add a new variable: 1) add it to the proper structure and make a copy in the
description structure 2) DescribeNetwork should initialize the description structure
3) GenerateNetwork should initialize the variable for each neuron/synapse


CAUTION: no saturation on external NMDA!!!
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <stdarg.h>
#include "rando2.h"


/*
Some stuff from Numerical Recipes
*/

#define IA 16807
#define IM 2147483647
#define AM (1.0/IM)
#define IQ 127773
#define IR 2836
#define NTAB 32
#define NDIV (1+(IM-1)/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)
//#define PI 3.141592654

long *idum;

void sran1(long *iseed)
{
  idum = iseed;
}

float ran1()
{
	int j;
	long k;
	static long iy=0;
	static long iv[NTAB];
	float temp;

	if (*idum <= 0 || !iy) {
		if (-(*idum) < 1) *idum=1;
		else *idum = -(*idum);
		for (j=NTAB+7;j>=0;j--) {
			k=(*idum)/IQ;
			*idum=IA*(*idum-k*IQ)-IR*k;
			if (*idum < 0) *idum += IM;
			if (j < NTAB) iv[j] = *idum;
		}
		iy=iv[0];
	}
	k=(*idum)/IQ;
	*idum=IA*(*idum-k*IQ)-IR*k;
	if (*idum < 0) *idum += IM;
	j=iy/NDIV;
	iy=iv[j];
	iv[j] = *idum;
	if ((temp=AM*iy) > RNMX) return RNMX;
	else return temp;
}


float gasdev()
{
	float ran1();
	static int iset=0;
	static float gset;
	float fac,rsq,v1,v2;

	if  (iset == 0) {
		do {
			v1=2.0*ran1()-1.0;
			v2=2.0*ran1()-1.0;
			rsq=v1*v1+v2*v2;
		} while (rsq >= 1.0 || rsq == 0.0);
		fac=sqrt(-2.0*log(rsq)/rsq);
		gset=v1*fac;
		iset=1;
		return v2*fac;
	} else {
		iset=0;
		return gset;
	}
}


//Below is the real stuff


#define MAXTN 10000 // max number of target neurons

// types of receptors

#define MAXRECEPTORS 3

#define ACh 0
#define GABAA 1
#define GABAB 2

// =============================================================
// Structures for the simulation

// --------------------- SYNAPSE structure ---------------------

typedef struct {
  int TargetPop;
  int TargetNeuron;
  int TargetAxon;
  float Efficacy; // change in the conductance due to a single spike (nS) (hence it is always positive, IPSP come out from rev pots)
  float LastConductance; // synaptic conductance when the last spike was emitted (useful for NMDA saturation), -1 if not NMDA (and hence no saturation)
  int TargetReceptor; // 0=Ach, 1=GABAA, 2=GABAB
  //Parameters and variables for the short term depression
  float D;                 // The fraction of available vesicle.
  float pv;                // Each spike reduce D by the factor pv.
  float tauD;             // The time constant for the speed of D recovery.

  //Parameters and variables for the short term facilitation
  float F;                 // The fraction of available vesicle.
  float Fp;                // Each spike increase F by the factor Fp.
  float tauF;             // The time constant for the speed of F decrease.

  // V. 7e9 for pre-synaptic inhibition
  float Ca_pre;   //Relative Ca level (max=1) at the presynaptic terminal
  float last_pre_update[MAXRECEPTORS];  // Time of the last update of Ca_pre
  float St_pre[MAXRECEPTORS];  // Gating variable on the pre-synaptic terminal
  float Ca_pre_tau[MAXRECEPTORS];        // Recovery time constant of Ca_pre (copied from the presynaptic time constant)

} Synapse;

// --------------------- NEURON structure -----------------------


typedef struct {

  // outputs

  int Naxonals;    // axonal tree (total number of synapses on the axon)
  Synapse *Axonals;

  // inputs

  int Nreceptors;
  float Tau[MAXRECEPTORS]; // tau for each conductance
  float LS[MAXRECEPTORS]; // total conductance for internal inputs (spikes generated by the network)
  float RevPots[MAXRECEPTORS]; // reversal potentials
  int   MgFlag[MAXRECEPTORS]; // 1 is magnesium block is active, 0 otherwise

  // external gaussian inputs for each receptor (to be added to S at each time step)
  float ExtS[MAXRECEPTORS]; // nS

  float ExtMuS[MAXRECEPTORS]; // nS/s
  float ExtSigmaS[MAXRECEPTORS]; // nS/s^2

  // here I should consider also the external inputs!!!

  // dynamic variables

  float V;
  int RefrState; // Refractory state counter
  float C; // capacitance in nF
  float Taum; // membrane time constant
  float RestPot; // Resting potential
  float ResetPot; // reset potential
  float Threshold; // threhsold
  float RefractT;    // Refractory period (ms)


  // Spike-frequency adaptation.
  // The properties of calcium-activated potassium current I

  float Ca;  //Concentration of Ca ions
  float Vk;  // resting potential
  float alpha_ca; // Amount of increment of [Ca] with each spike discharge. (muM)
  float tau_ca; // time constant
  float g_ahp; // efficacy

  // the times of the last two spikes
  float TimeLastSpike; // time of the last emitted spike (for NMDA saturation)
  float PTimeLastSpike;



  //Anomalous delayed rectifier (ADR). Simulating the up and down state in striatal neurons.
  float g_adr_max;  //Maximun value of the g
  float Vadr_h; //Potential for g_adr=0.5g_adr_max
  float Vadr_s; //Slop of g_adr at Vadr_h, defining how sharp the shape of g_ard is.
  float ADRRevPot; //reverse potential for ADR.
  float g_k_max;  // Maximun outward rectifying current
  float Vk_h; //Potential for g_k=0.5g_k_max
  float Vk_s; //Slop of g_k at Vk_h, defining how sharp the shape of g_k is.
  float tau_k_max; //maximum value of tau_k
  float n_k;  // gating variable for outward rectifying K channel;


} Neuron;



// ------------------- Population structure ---------------------
#define MAXDELAYINDT 50
#define MAXSPIKESINDT 2000

typedef struct {
  char Label[100];

  int Ncells;
  Neuron *Cell;

  // Table of spikes emitted by the neurons of this population
  int CTableofSpikes; // pointer where we are (time t)
  int DTableofSpikes; // pointer to the t-delay
  int NTableofSpikes[MAXDELAYINDT];
  int TableofSpikes[MAXDELAYINDT][MAXSPIKESINDT];

} Population;

#define MAXP 300
#define MAXEXTMEM 20

int Npop;
Population Pop[MAXP];

// global variables for the simulation

float dt=0.1;      // in ms
float Time;    // absolute time from the beginning of the trial
int delayindt; // transmission delay in dt
int flagverbose=1; // flag to activate verbose messages
int FlagSaveAllSpikes=1;
float ext_freq; // external noise
// =================================================================
// Descriptors to generate the network structure

typedef struct {
  float Connectivity; // mean fraction of randomly connected target neurons
  float TargetReceptor; // 0=Ach, ...
  float MeanEfficacy; // mean efficacy (for initialization)
  float EfficacySD; // Standard deviation of the efficacy distribution.
  float pv;  // Each spike reduce the fraction D of available vesicle by the factor pv.
  float tauD; // The time constant for the speed of D recovery.
  float Fp;   // Each spike increase F by the factor Fp.
  float tauF; // The time constant for the speed of F decrease.
//  float Ca_pre_tau;  // Time time constant for the pre-synaptic inhibition
} SynPopDescr;

typedef struct {
  char Label[100];
  int Ncells;

  int NTargetPops;
  int TargetPops[MAXP];
  int TargetAxons[MAXP];

  SynPopDescr SynP[MAXP];

  int Nreceptors;
  char  ReceptorLabel[MAXRECEPTORS][100]; // label for the receptor (needed for the compiler)
  float Tau[MAXRECEPTORS]; // tau for each conductance
  float RevPots[MAXRECEPTORS]; // reversal potentials
  int   MgFlag[MAXRECEPTORS]; // magnesium block flag

  // external input (externally defined)
  float MeanExtCon[MAXRECEPTORS]; // mean total number of external connections
  float MeanExtEff[MAXRECEPTORS]; // external mean efficacy
  float ExtEffSD[MAXRECEPTORS]; // Deviation of distribution of external efficacy
  float FreqExt[MAXRECEPTORS]; // external frequency in Hz
 float FreqExtTau[MAXRECEPTORS]; // time constant of the decaly from old FreqExt to new FreqExt
  float FreqExtSD[MAXRECEPTORS]; // external frequency in Hz
  float FreqExtMem[MAXRECEPTORS]; // memory in the external noise
  float FreqExtDecay[MAXRECEPTORS]; //decay factor of the external noise
  float FreqExtNorm[MAXRECEPTORS]; //Normalization factor for the memory of the external noise

  // external input (statistics internally calculated)
  float MeanExtMuS[MAXRECEPTORS]; // statistics of the external input nS/s
  float MeanExtSigmaS[MAXRECEPTORS];
  // dynamic variables

  float C; // capacitance
  float Taum; // membrane time constant
  float RestPot; // Resting potential
  float ResetPot; // reset potential
  float Threshold; // threhsold
  float RefractT; // Refractory period (ms)

  float Vk;  // resting potential
  float alpha_ca; // Amount of increment of [Ca] with each spike discharge. (muM)
  float tau_ca; // time constant
  float g_ahp; // efficacy

  //Anomalous delayed rectifier (ADR)
  float g_adr_max;  //Maximun value of the g
  float Vadr_h; //Potential for g_adr=0.5g_adr_max
  float Vadr_s; //Slop of g_adr at Vadr_h, defining how sharp the shape of g_ard is.
  float ADRRevPot; //Reverse potential for ADR
  float g_k_max;  // Maximun outward rectifying current
  float Vk_h; //Potential for g_k=0.5g_k_max
  float Vk_s; //Slop of g_k at Vk_h, defining how sharp the shape of g_k is.
  float tau_k_max; //maximum tau for outward rectifying k current


} PopDescr;

PopDescr PopD[MAXP];

// ===============================================================
// Protocol descriptors

#define MAXEVENTS 200

#define ENDOFTRIAL 1
#define CHANGEEXTFREQ 2
#define CHANGEMEANEFF 3
#define CHANGEEXTFREQSD 4

typedef struct {
  int Type;
  float ETime;
  int PopNumber;
  int TargetPopNumber;
  int ReceptorNumber;
  float FreqExt;
  float FreqExtDecay;
  float FreqExtSD;
  float FreqExtsd;
  float MeanEff;
  char Label[100];
} EventDescr;

int NEvents=0;
int CEvent; // current event
float NextEventTime=0.;
float TrialDuration=1000.; // in ms
EventDescr Events[MAXEVENTS];

// Multitrial version:

int CurrentTrial;
int NumberofTrials;


// ===============================================================================
// AUXILIARY SUBROUTINES
// ===============================================================================

/* --------------------------------------------------------------------------------
   Report: prints on screen and log file (dev_oss filename device)
   -------------------------------------------------------------------------------- */

void report(char *fmt,...)
{
  static FILE *devlog;
  static int initflag=1;
  va_list ap;
  char *s,fmtemp[20];
  int ival,fmtemp_i;
  float fval;
  char *cval;

  if(initflag) {
    devlog=fopen("simu.log","w");
    if(devlog==NULL) return;
    initflag=0;
  }

  va_start(ap,fmt);
  for(s=fmt; *s; s++) {

    if(*s!='%') { if(flagverbose) printf("%c",*s); fprintf(devlog,"%c",*s); continue; }

    fmtemp_i=0; while (*s && !(*s=='s' || *s=='d' || *s=='f')) { fmtemp[fmtemp_i]=*s; s++; fmtemp_i++; }
    if(*s=='d') {  ival=va_arg(ap,int); fmtemp[fmtemp_i]='d'; fmtemp[fmtemp_i+1]=0;
                   if(flagverbose) printf(fmtemp,ival); fprintf(devlog,fmtemp,ival);  }
    if(*s=='f') {  fval=va_arg(ap,double); fmtemp[fmtemp_i]='f'; fmtemp[fmtemp_i+1]=0;
                   if(flagverbose) printf(fmtemp,fval); fprintf(devlog,fmtemp,fval);  }
    if(*s=='s') {  cval=va_arg(ap,char *); fmtemp[fmtemp_i]='s'; fmtemp[fmtemp_i+1]=0;
                   if(flagverbose) printf(fmtemp,cval); fprintf(devlog,fmtemp,cval);  }
  } /* end for */
  va_end(ap);
}



// ===============================================================
//
// INPUT
// ==============================================================
// auxiliary subroutines for parsing the descriptor file
// returns the code of the population with name s (-1 in case of error)

int PopulationCode(char *s)
{
  int p;


  for(p=0;p<Npop;p++)
    {

      if(strcmp(PopD[p].Label,s)==0) {
	return p;
      }
    }
  return -1;
}


// returns the code of the receptor with name s for pop p (-1 in case of error)

int ReceptorCode(char *s,int p)
{
  int r;

  for(r=0;r<PopD[p].Nreceptors;r++)
    {
      if(strcmp(PopD[p].ReceptorLabel[r],s)==0) {
	return r;
      }
    }
  return -1;
}

// =======================================================================
// DescribeNetwork()
// Initializes the descriptors of the network by parsing cl_network.conf
// =======================================================================

#define EXC 0
#define INH 1


int DescribeNetwork()
{

  FILE *devconf;
  char buf[1000],*s,*es;
  int currentpopflag=0;
  int currentreceptorflag=0;
  int line,auxi;
  int currentpop,currentreceptor,currenttarget,currenttargetaxon;
  float aux;
  float std_p, std_tau;
  int flag_d;

  // FIRST PASSAGE
  // -------------------------------------------------------------------
  // the parser has to go over the file twice: the first time reads all
  // the labels and initializes the number of populations and the number
  // of receptors
  report("Parsing network configuration... first passage\n");

  /*  strncpy=(network_conf,prefix,strlen(prefix));
  strcpy=(network_conf+strlen(prefix),"network.conf");
  devconf=fopen(network_conf,"r");*/
  devconf=fopen("network.conf","r");
  if(devconf==NULL) { printf("ERROR:  Unable to read configuration file\n"); return 0;}

  Npop=0; line=-1;

  while(fgets(buf,1000,devconf)!=NULL)
    {
      line++;
      s=buf;
      // trim \n at the end
      es=buf; while(*es && *es!='\n') es++; *es=0;

      while(*s==' ' || *s=='\t') s++; // skip blanks
      if(*s==0) continue; // empty line
      if(*s=='%') continue; // remark

      // commands for defining a new population
      if(strncmp(s,"NeuralPopulation:",17)==0)
	{
	  currentpopflag=1;
	  s+=17; while(*s==' ') s++;
	  strcpy(PopD[Npop].Label,s);
	  report("Population: %s\n",PopD[Npop].Label);
	  PopD[Npop].Nreceptors=0;
	  continue;
	}

      if(strncmp(s,"EndNeuralPopulation",19)==0)
	{
	  currentpopflag=0;
	  Npop++;
	  continue;
	}

      // command for defining a receptor
      if(strncmp(s,"Receptor:",9)==0) {
	currentreceptorflag=1;
	s+=9; while(*s==' ') s++;
	strcpy(PopD[Npop].ReceptorLabel[PopD[Npop].Nreceptors],s);
	report("Receptor %d: %s\n",PopD[Npop].Nreceptors,PopD[Npop].ReceptorLabel[PopD[Npop].Nreceptors]);
      }

      if(strncmp(s,"EndReceptor",11)==0) {
	currentreceptorflag=0;
	PopD[Npop].Nreceptors++;
      }
    }

  fclose(devconf);

  // Second passage: now all the parameters and the target populations are parsed
  // ----------------------------------------------------------------------------

  report("Parsing network configuration... second passage\n");

  //  devconf=fopen(network_conf,"r");
  devconf=fopen("network.conf","r");
  if(devconf==NULL) { printf("ERROR:  Unable to read configuration file\n"); return 0;}

  line=-1;

  while(fgets(buf,1000,devconf)!=NULL)
    {
      line++;
      s=buf;

      // trim \n at the end
      es=buf; while(*es && *es!='\n') es++; *es=0;

      while(*s==' ' || *s=='\t') s++; // skip blanks
      if(*s==0) continue; // empty line
      if(*s=='%') continue; // remark

      // population commands
      if(strncmp(s,"NeuralPopulation:",17)==0)
	{
	  currentpopflag=1;
	  s+=17; while(*s==' ') s++;
	  currentpop=PopulationCode(s);
	  if(currentpop==-1) { printf("Unknown population [%s]: line %d\n",s,line); return -1;}
	  PopD[currentpop].NTargetPops=0;
	  report("-----------------------------------\n    Population: %s\n-----------------------------------\n",PopD[currentpop].Label);

	  //Initialize some population parameters
         PopD[currentpop].RefractT=2.0;
          PopD[currentpop].g_adr_max=0;
          PopD[currentpop].Vadr_h=-100;
          PopD[currentpop].Vadr_s=10;
          PopD[currentpop].ADRRevPot=-90;
          PopD[currentpop].g_k_max=0;
          PopD[currentpop].Vk_h=-34;
	  PopD[currentpop].Vk_s=6.5;
          PopD[currentpop].tau_k_max=8;
	  continue;
	}

      if(strncmp(s,"EndNeuralPopulation",19)==0)
	{
	  report("EndPopulation\n");
	  currentpopflag=0;
	  continue;
	}

      // paramters for the population


      if(strncmp(s,"N=",2)==0 && currentpopflag) {
	PopD[currentpop].Ncells=atoi(s+2); report("  N=%d\n",PopD[currentpop].Ncells); continue;
      }
      if(strncmp(s,"C=",2)==0 && currentpopflag) {
	PopD[currentpop].C=atof(s+2);report("  C=%f nF\n",(double)PopD[currentpop].C); continue;
      }
      if(strncmp(s,"Taum=",5)==0 && currentpopflag) {
	PopD[currentpop].Taum=atof(s+5);report("  Membrane time constant=%f ms\n",(double)PopD[currentpop].Taum); continue;
      }
      if(strncmp(s,"RestPot=",8)==0 && currentpopflag) {
	PopD[currentpop].RestPot=atof(s+8);report("  Resting potential=%f mV\n",(double)PopD[currentpop].RestPot); continue;
      }
      if(strncmp(s,"ResetPot=",9)==0 && currentpopflag) {
	PopD[currentpop].ResetPot=atof(s+9);report("  Reset potential=%f mV\n",(double)PopD[currentpop].ResetPot); continue;
      }
      if(strncmp(s,"Threshold=",10)==0 && currentpopflag) {
	PopD[currentpop].Threshold=atof(s+10);report("  Threshold =%f mV\n",(double)PopD[currentpop].Threshold); continue;
      }

      if(strncmp(s,"RestPot_ca=",11)==0 && currentpopflag) {
	PopD[currentpop].Vk=atof(s+11);report("  RestPot_ca =%f mV\n",(double)PopD[currentpop].Vk); continue;
      }

      if(strncmp(s,"RefractT=",9)==0 && currentpopflag) {
	PopD[currentpop].RefractT=atof(s+9);report("  Refractory period =%f mV\n",(double)PopD[currentpop].RefractT); continue;
      }

      if(strncmp(s,"Alpha_ca=",9)==0 && currentpopflag) {
	PopD[currentpop].alpha_ca=atof(s+9);report("  Alpha_ca =%f mV\n",(double)PopD[currentpop].alpha_ca); continue;
      }

      if(strncmp(s,"Tau_ca=",7)==0 && currentpopflag) {
	PopD[currentpop].tau_ca=atof(s+7);report("  Tau_ca =%f mV\n",(double)PopD[currentpop].tau_ca); continue;
      }

      if(strncmp(s,"Eff_ca=",7)==0 && currentpopflag) {
	PopD[currentpop].g_ahp=atof(s+7);report("  Eff_ca =%f mV\n",(double)PopD[currentpop].g_ahp); continue;
      }

      if(strncmp(s,"g_ADR_Max=",10)==0 && currentpopflag) {
	PopD[currentpop].g_adr_max=atof(s+10);report("  g_adr_max=%f mV\n",(double)PopD[currentpop].g_adr_max); continue;
      }

      if(strncmp(s,"V_ADR_h=",8)==0 && currentpopflag) {
	PopD[currentpop].Vadr_h=atof(s+8);report("  V_ADR_h=%f mV\n",(double)PopD[currentpop].Vadr_h); continue;
      }

      if(strncmp(s,"V_ADR_s=",8)==0 && currentpopflag) {
	PopD[currentpop].Vadr_s=atof(s+8);report("  V_ADR_s=%f mV\n",(double)PopD[currentpop].Vadr_h); continue;
      }

      if(strncmp(s,"ADRRevPot=",10)==0 && currentpopflag) {
	PopD[currentpop].ADRRevPot=atof(s+10);report("  ADRRevPot=%f mV\n",(double)PopD[currentpop].ADRRevPot); continue;
      }

      if(strncmp(s,"g_k_Max=",8)==0 && currentpopflag) {
	PopD[currentpop].g_k_max=atof(s+8);report("  g_k_max=%f mV\n",(double)PopD[currentpop].g_k_max); continue;
      }

      if(strncmp(s,"V_k_h=",6)==0 && currentpopflag) {
	PopD[currentpop].Vk_h=atof(s+6);report("  V_k_h=%f mV\n",(double)PopD[currentpop].Vk_h); continue;
      }

      if(strncmp(s,"V_k_s=",6)==0 && currentpopflag) {
      	PopD[currentpop].Vk_s=atof(s+6);report("  V_k_s=%f mV\n",(double)PopD[currentpop].Vk_h); continue;
         }
      if(strncmp(s,"tau_k_max=",10)==0 && currentpopflag) {
	PopD[currentpop].tau_k_max=atof(s+10);report("  tau_k_max=%f mV\n",(double)PopD[currentpop].tau_k_max); continue;
      }

      // receptor paramters
      if(strncmp(s,"Receptor:",9)==0 && currentpopflag) {
	currentreceptorflag=1;
	s+=9; while(*s==' ') s++;
	currentreceptor=ReceptorCode(s,currentpop);
	if(currentreceptor==-1) { printf("ERROR: Unknown receptor type\n"); return -1; }
	if(strncmp(s,"GABAB",4)==0) { // deactivate magnesium block for GABAB
	  PopD[currentpop].MgFlag[currentreceptor]=0;
	}
	else PopD[currentpop].MgFlag[currentreceptor]=0;
	report("Receptor %d: %s (Mg block: %d)\n",currentreceptor,PopD[currentpop].ReceptorLabel[currentreceptor],PopD[currentpop].MgFlag[currentreceptor]);

        PopD[currentpop].ExtEffSD[currentreceptor]=0;
        PopD[currentpop].FreqExtSD[currentreceptor]=0;
	continue;
      }

      if(strncmp(s,"EndReceptor",11)==0) {
	currentreceptorflag=0;
	continue;
      }

      if(strncmp(s,"Tau=",4)==0 && currentpopflag && currentreceptorflag) {
	PopD[currentpop].Tau[currentreceptor]=atof(s+4);report("  Tau=%f ms\n",(double)PopD[currentpop].Tau[currentreceptor]); continue;
      }
      if(strncmp(s,"RevPot=",7)==0 && currentpopflag && currentreceptorflag) {
	PopD[currentpop].RevPots[currentreceptor]=atof(s+7);report("  Reversal potential=%f mV\n",(double)PopD[currentpop].RevPots[currentreceptor]); continue;
      }
      if(strncmp(s,"FreqExt=",8)==0 && currentpopflag && currentreceptorflag) {
	PopD[currentpop].FreqExt[currentreceptor]=atof(s+8);report("  Ext frequency=%f Hz\n",(double)PopD[currentpop].FreqExt[currentreceptor]); continue;
      }
      if(strncmp(s,"FreqExtSD=",10)==0 && currentpopflag && currentreceptorflag) {
	PopD[currentpop].FreqExtSD[currentreceptor]=atof(s+10);report("  Ext frequency SD=%f Hz\n",(double)PopD[currentpop].FreqExtSD[currentreceptor]); continue;
      }
      if(strncmp(s,"MeanExtEff=",11)==0 && currentpopflag && currentreceptorflag) {
	PopD[currentpop].MeanExtEff[currentreceptor]=atof(s+11);report("  Ext efficacy=%f nS\n",(double)PopD[currentpop].MeanExtEff[currentreceptor]); continue;}
      if(strncmp(s,"ExtEffSD=",9)==0 && currentpopflag && currentreceptorflag) {
	PopD[currentpop].ExtEffSD[currentreceptor]=atof(s+9);report("  Ext efficacy SD=%f nS\n",(double)PopD[currentpop].ExtEffSD[currentreceptor]); continue;
      }
      if(strncmp(s,"MeanExtCon=",11)==0 && currentpopflag && currentreceptorflag) {
	PopD[currentpop].MeanExtCon[currentreceptor]=atof(s+11);report("  Ext connections=%f\n",(double)PopD[currentpop].MeanExtCon[currentreceptor]); continue;
      }

      // target populations

      if(strncmp(s,"TargetPopulation:",17)==0 && currentpopflag)
	{
	  s+=17; while(*s==' ') s++;
	  currenttarget=PopulationCode(s);
	  if(currenttarget==-1) { printf("Unknown target population: line %d\n",line); return -1;}
	  PopD[currentpop].TargetPops[PopD[currentpop].NTargetPops]=currenttarget;

	  report("Synapses targeting population: %s (%d)\n ",PopD[currenttarget].Label,currenttarget);

          PopD[currentpop].TargetAxons[PopD[currentpop].NTargetPops]=-1;
          PopD[currentpop].SynP[PopD[currentpop].NTargetPops].pv=0;
          PopD[currentpop].SynP[PopD[currentpop].NTargetPops].tauD=1000;
          PopD[currentpop].SynP[PopD[currentpop].NTargetPops].Fp=0;
          PopD[currentpop].SynP[PopD[currentpop].NTargetPops].tauF=5000;
          PopD[currentpop].SynP[PopD[currentpop].NTargetPops].EfficacySD=0;
//          PopD[currentpop].SynP[PopD[currentpop].NTargetPops].Ca_pre_tau=dt;
	  continue;
	}

      if(strncmp(s,"EndTargetPopulation",20)==0)
	{
	  PopD[currentpop].NTargetPops++;
	  continue;
	}

      if(strncmp(s,"TargetAxon=",11)==0 && currentpopflag)
	{
	  s+=11; while(*s==' ') s++;
	  currenttargetaxon=PopulationCode(s);
	  if(currenttarget==-1) { printf("Unknown target population: line %d\n",line); return -1;}
	  PopD[currentpop].TargetAxons[PopD[currentpop].NTargetPops]=currenttargetaxon;
          report("  TargetAxon=%s\n",s);
	}
      if(strncmp(s,"Connectivity=",13)==0 && currentpopflag)
	{
	  aux=atof(s+13);
	  PopD[currentpop].SynP[PopD[currentpop].NTargetPops].Connectivity=aux;report("  Connectivity=%f\n",(double)aux);
	}
      if(strncmp(s,"TargetReceptor=",15)==0 && currentpopflag)
	{
	  s+=15; while(*s==' ' || *s=='\t') s++;
	  auxi=ReceptorCode(s,currenttarget);
	  if(auxi==-1) { printf("Unknown target receptor, line %d\n",line); return -1; }
	  PopD[currentpop].SynP[PopD[currentpop].NTargetPops].TargetReceptor=auxi;
	  report("  Target receptor code=%d\n",auxi);
	}

      if(strncmp(s,"MeanEff=",8)==0 && currentpopflag)
	{
	  aux=atof(s+8);
	  PopD[currentpop].SynP[PopD[currentpop].NTargetPops].MeanEfficacy=aux;report("  Mean efficacy=%f\n",(double)aux);
	}
      if(strncmp(s,"EffSD=",6)==0 && currentpopflag)
	{
	  aux=atof(s+6);
	  PopD[currentpop].SynP[PopD[currentpop].NTargetPops].EfficacySD=aux;report("  Efficacy SD=%f\n",(double)aux);
	}

      if(strncmp(s,"STDepressionP=",14)==0 && currentpopflag) {
        aux=atof(s+14);
	PopD[currentpop].SynP[PopD[currentpop].NTargetPops].pv=aux;
        report("  STDepressionP=%f mV\n",(double)aux);
      }

      if(strncmp(s,"STDepressionTau=",16)==0 && currentpopflag) {
        aux=atof(s+16);
	PopD[currentpop].SynP[PopD[currentpop].NTargetPops].tauD=aux;
        report("  STDepressionTau=%f mV\n",(double)aux);
      }
      if(strncmp(s,"STFacilitationP=",16)==0 && currentpopflag) {
        aux=atof(s+16);
	PopD[currentpop].SynP[PopD[currentpop].NTargetPops].Fp=aux;
        report("  STFacilitationP=%f mV\n",(double)aux);
      }

      if(strncmp(s,"STFacilitationTau=",18)==0 && currentpopflag) {
        aux=atof(s+18);
	PopD[currentpop].SynP[PopD[currentpop].NTargetPops].tauF=aux;
        report("  STFacilitationTau=%f mV\n",(double)aux);
      }
    } // end while


  fclose(devconf);

 	return 1;

}


// -------------------------------------------------------------------------------------
// Parse protocol
// -------------------------------------------------------------------------------------

int ParseProtocol()
{
  FILE *devprot;
  char buf[1000],*s,*es;
  int eventflag=0;
  int line,auxi,currentpop;
  float aux;

  report("-------------------------------------------------\nParsing protocol...\n");

  /*  strncpy=(network_pro,prefix,strlen(prefix));
  strcpy=(network_pro+strlen(prefix),"network.pro");
  devconf=fopen(network_pro,"r");*/
  devprot=fopen("network.pro","r");
  if(devprot==NULL) { printf("ERROR:  Unable to read protocol file\n"); return 0;}

  line=-1; NEvents=0;

  while(fgets(buf,1000,devprot)!=NULL)
    {
      line++;
      s=buf;
      // trim \n at the end
      es=buf; while(*es && *es!='\n') es++; *es=0;

      while(*s==' ' || *s=='\t') s++; // skip blanks
      if(*s==0) continue; // empty line
      if(*s=='%') continue; // remark



      // commands for defining a new event
      if(strncmp(s,"EventTime",9)==0)
	{
	  eventflag=1;
	  s+=9;
	  Events[NEvents].ETime=atof(s);
	  if(Events[NEvents].ETime<0.)
	    {
	      printf("ERROR: Invalid event time, line %d\n",line);
	      return -1;
	    }
          Events[NEvents].FreqExtsd=0;    // initialize FreqExtsd
	  report("Event %d: time %f ms\n",NEvents,Events[NEvents].ETime);
	  continue;
	}

      if(strncmp(s,"EndEvent",8)==0)
	{
	  eventflag=0;
	  NEvents++;
	  continue;
	}

      if(strncmp(s,"Type=",5)==0) {
	s+=5; while(*s==' ' || *s=='\t') s++;
	if(strncmp(s,"ChangeExtFreq",13)==0) {
	  Events[NEvents].Type=CHANGEEXTFREQ;
         Events[NEvents].FreqExtDecay=dt;
	}
	if(strncmp(s,"ChangeExtFreqSD",15)==0) {
	  Events[NEvents].Type=CHANGEEXTFREQSD;
	}
	if(strncmp(s,"ChangeMeanEff",13)==0) {
	  Events[NEvents].Type=CHANGEMEANEFF;
	}
	if(strncmp(s,"EndTrial",8)==0) {
	  Events[NEvents].Type=ENDOFTRIAL;
	  TrialDuration=Events[NEvents].ETime;
	}
	continue;
      }

      if(strncmp(s,"FreqExt=",8)==0) {
	s+=8;
	Events[NEvents].FreqExt=atof(s);
	report("  New external frequency: %f Hz\n",Events[NEvents].FreqExt);
	continue;
      }
      if(strncmp(s,"FreqExtDecay=",13)==0) {
	s+=13;
	Events[NEvents].FreqExtDecay=atof(s);
	report("  Time constant for reaching the new external frequency: %f Hz\n",Events[NEvents].FreqExtDecay);
	continue;
      }
      if(strncmp(s,"FreqExtSD=",10)==0) {
	s+=10;
	Events[NEvents].FreqExtSD=atof(s);
	report("  New external frequency SD: %f Hz\n",Events[NEvents].FreqExtSD);
	continue;
      }
      if(strncmp(s,"FreqExtsd=",10)==0) {
	s+=10;
	Events[NEvents].FreqExtsd=atof(s);
	report("  New external frequency offset: %f Hz\n",Events[NEvents].FreqExtsd);
	continue;
      }
      if(strncmp(s,"MeanEff=",8)==0) {
	s+=8;
	Events[NEvents].MeanEff=atof(s);
	report("  New Mean Eff: %f Hz\n",Events[NEvents].MeanEff);
	continue;
      }

      if(strncmp(s,"Label=",6)==0) {
	s+=6;
	strcpy(Events[NEvents].Label,s);
	report("   Label: [%s]\n",Events[NEvents].Label);
	continue;

      }

      if(strncmp(s,"Population:",11)==0 && eventflag) {
	s+=11;	while(*s==' ') s++;
	currentpop=PopulationCode(s);
	if(currentpop==-1) { printf("ERROR: Unknown population: line %d\n",line); return -1;}
	Events[NEvents].PopNumber=currentpop;
	report("  Population code: %d\n",currentpop);
	continue;
      }

      if(strncmp(s,"TargetPopulation:",17)==0 && eventflag) {
	s+=17;	while(*s==' ') s++;
	currentpop=PopulationCode(s);
	if(currentpop==-1) { printf("ERROR: Unknown population: line %d\n",line); return -1;}
	Events[NEvents].TargetPopNumber=currentpop;
	report("  TargetPopulation code: %d\n",currentpop);
	continue;
      }

      if(strncmp(s,"Receptor:",9)==0 && eventflag)
	{
	  s+=9; while(*s==' ' || *s=='\t') s++;
	  auxi=ReceptorCode(s,currentpop);
	  if(auxi==-1) { printf("ERROR: Unknown receptor, line %d\n",line); return -1; }
	  Events[NEvents].ReceptorNumber=auxi;
	  report("  Receptor code:%d\n",auxi);

	}

    } // end while

  // sort events!... (for now I rely on the fact that events are sorted)

  return 1;

}


// ======================================================================================
// GenerateNetwork()
//
// Generates all the structures which are needed to run the simulation
// on the basis of PopD structures. It also allocates the memory for the
// network. PopD is initialized in DescribeNetwork.
// ======================================================================================

int GenerateNetwork()
{

  int p,i,j,r,tp,tpi,tn,k,rd,na,tpa,flag_a;
  int ni[MAXTN],tpni[MAXTN],trni[MAXTN], tpax[MAXTN];
  float eff[MAXTN],ts[MAXTN],timets[MAXTN],D[MAXTN],pv[MAXTN],tauD[MAXTN],F[MAXTN],Fp[MAXTN],tauF[MAXTN];
  float efficacy;

  if(flagverbose) {
    report("\nGenerating network...\n");
  }


  delayindt=1;

  // loop over all the populations
  for(p=0;p<Npop;p++)
    {
      strcpy(Pop[p].Label,PopD[p].Label);
      if(flagverbose) {
	report("Population %2d: %s\n",p,Pop[p].Label);
      }

      for(rd=0;rd<PopD[p].Nreceptors;rd++)
	{
	  // external input: compute first the asymptoic mu and sigma (m,s in nphys notation) of the conductances
	  //          do{efficacy=(gasdev()*PopD[p].ExtEffSD[rd])+PopD[p].MeanExtEff[rd];}while(efficacy<0);
          efficacy=PopD[p].MeanExtEff[rd];
	  PopD[p].MeanExtMuS[rd]=PopD[p].FreqExt[rd]*.001*efficacy*PopD[p].MeanExtCon[rd]*PopD[p].Tau[rd];
	  PopD[p].MeanExtSigmaS[rd]=sqrt(PopD[p].Tau[rd]*.5*PopD[p].FreqExt[rd]*.001*efficacy*efficacy*PopD[p].MeanExtCon[rd]);
          PopD[p].FreqExtTau[rd]=dt;
          PopD[p].FreqExtDecay[rd]=0.999;
          PopD[p].FreqExtNorm[rd]=PopD[p].FreqExtDecay[rd]/(1-PopD[p].FreqExtDecay[rd]);
          PopD[p].FreqExtMem[rd]=PopD[p].MeanExtEff[rd]*PopD[p].FreqExtNorm[rd];
	  report("Pop %d Conductance: %d m=%f nS s=%f nS\n",p,rd,PopD[p].MeanExtMuS[rd],PopD[p].MeanExtSigmaS[rd]);
	}


      // allocate memory for all neurons
      Pop[p].Ncells=PopD[p].Ncells;
      if(flagverbose) { report("  - Number of cells: %d\n",Pop[p].Ncells); }
      Pop[p].Cell=(Neuron *) calloc(Pop[p].Ncells,sizeof(Neuron));

      // init auxiliary variables like the tables of spikes
      if(delayindt>MAXDELAYINDT) { printf("ERROR: delay too long. Increase MAXDELAYINDT and recompile\n"); return -1;}
      Pop[p].CTableofSpikes=delayindt;
      Pop[p].DTableofSpikes=0;
      for(k=0;k<MAXDELAYINDT;k++)
	{
	  Pop[p].NTableofSpikes[k]=0;
	}

      // generate neurons and their axonal trees
      for(i=0;i<Pop[p].Ncells;i++)
	{
	  // single neuron parameters
	  Pop[p].Cell[i].Taum=PopD[p].Taum;
	  Pop[p].Cell[i].ResetPot=PopD[p].ResetPot;
	  Pop[p].Cell[i].C=PopD[p].C;
	  Pop[p].Cell[i].RestPot=PopD[p].RestPot;
	  Pop[p].Cell[i].Threshold=PopD[p].Threshold;
	  Pop[p].Cell[i].RefractT=PopD[p].RefractT;
	  Pop[p].Cell[i].Vk=PopD[p].Vk;
	  Pop[p].Cell[i].alpha_ca=PopD[p].alpha_ca;
	  Pop[p].Cell[i].tau_ca=PopD[p].tau_ca;
	  Pop[p].Cell[i].g_ahp=PopD[p].g_ahp;
          Pop[p].Cell[i].g_adr_max=PopD[p].g_adr_max;
          Pop[p].Cell[i].Vadr_h=PopD[p].Vadr_h;
          Pop[p].Cell[i].Vadr_s=PopD[p].Vadr_s;
          Pop[p].Cell[i].ADRRevPot=PopD[p].ADRRevPot;
          Pop[p].Cell[i].g_k_max=PopD[p].g_k_max;
          Pop[p].Cell[i].Vk_h=PopD[p].Vk_h;
          Pop[p].Cell[i].Vk_s=PopD[p].Vk_s;
          Pop[p].Cell[i].tau_k_max=PopD[p].tau_k_max;
          Pop[p].Cell[i].n_k=0;


	  // receptors
	  Pop[p].Cell[i].Nreceptors=PopD[p].Nreceptors;
	  for(r=0;r<Pop[p].Cell[i].Nreceptors;r++)
	    {
	      // parameters (simply copy if the network is homoegeneous)
	      Pop[p].Cell[i].Tau[r]=PopD[p].Tau[r];
	      Pop[p].Cell[i].RevPots[r]=PopD[p].RevPots[r];
	      Pop[p].Cell[i].MgFlag[r]=PopD[p].MgFlag[r]; // magnesium block
	      // init
	      Pop[p].Cell[i].LS[r]=0.;

	      // external input
	      do{efficacy=(gasdev()*PopD[p].ExtEffSD[r])+PopD[p].MeanExtEff[r];}while(efficacy<0);
  	      Pop[p].Cell[i].ExtMuS[r]=PopD[p].FreqExt[r]*.001*efficacy*PopD[p].MeanExtCon[r]*PopD[p].Tau[r];
	      Pop[p].Cell[i].ExtSigmaS[r]=sqrt(PopD[p].Tau[r]*.5*PopD[p].FreqExt[r]*.001*efficacy*efficacy*PopD[p].MeanExtCon[r]);
	      //if(r==0)printf("%f\n",Pop[p].Cell[i].ExtMuS[r]);
	      Pop[p].Cell[i].ExtS[r]=0.;

	    } // end for r

	  // init
	  Pop[p].Cell[i].V=Pop[p].Cell[i].RestPot;
	  Pop[p].Cell[i].RefrState=0;
	  Pop[p].Cell[i].TimeLastSpike=-10000.; // time of the last emitted spike (for NMDA saturation)
	  Pop[p].Cell[i].PTimeLastSpike=-10000.; // time of spike emitted previous the last emitted spiek (for NMDA saturation)

	  // Axonal tree -------------------------------
	  // loop on target populations
	  Pop[p].Cell[i].Naxonals=0;

     	  for(tpi=0;tpi<PopD[p].NTargetPops;tpi++)
	    {
	      tp=PopD[p].TargetPops[tpi]; // target population (tpi=target pop index)
	      // loop on the neurons of the target population

	      // choose the target neurons/populations and put them in a temporary vector ni/tpni
	      na=Pop[p].Cell[i].Naxonals; // auxiliary variable to speed up

              // targeted axon (for pre-synaptic inhibition)
              tpa=PopD[p].TargetAxons[tpi];


              // Source and target populations have to have the same number of neurons to be
              // connected one-to-one
              if(PopD[tp].Ncells!=PopD[p].Ncells){printf("Error: Source and target populations do not have the same number of neurons. Cannot make one-to-one connection. Exit..."); exit(0);}
              //Neuron cannot connect to itself
              if(tp==p){printf("Error: Neuron cannot connect to itself. Exit..."); exit(0);}

              // To make one-to-one connection, connectivity should be equal to -1
              if(PopD[p].SynP[tpi].Connectivity!=-1)
                {printf("Error: To make one-to-one connection, the value of connectivity should be set to -1. Exit..."); exit(0);}


//	      for(tn=0;tn<PopD[tp].Ncells;tn++)
//		{


//		  if(drand49()<PopD[p].SynP[tpi].Connectivity) { // the connection exists
		    ni[na]=i;   // target neuron
		    tpni[na]=tp; // target population
                    tpax[na]=tpa;
		    trni[na]=PopD[p].SynP[tpi].TargetReceptor; // target receptor
                    do{efficacy=gasdev()*PopD[p].SynP[tpi].EfficacySD+PopD[p].SynP[tpi].MeanEfficacy;}while(efficacy<0);
		    eff[na]= efficacy;// efficacy
		    // printf("%f\n",efficacy);
                    D[na]=1;
                    pv[na]=PopD[p].SynP[tpi].pv;
                    tauD[na]=PopD[p].SynP[tpi].tauD;
                    Fp[na]=PopD[p].SynP[tpi].Fp;
                    if(Fp[na]==0)
                       F[na]=1;
                    else
                       F[na]=0;
                    tauF[na]=PopD[p].SynP[tpi].tauF;

		    if(strncmp(PopD[tp].ReceptorLabel[trni[na]],"NMDA",4)==0) {
		      ts[na]=0.; // initial LastConductance (for NMDA saturation)
		    }
		    else ts[na]=-1.; // not an NMDA type (so, no saturation)
		    na++;
		    Pop[p].Cell[i].Naxonals++;
//		  }
//		} // end for tn
	    } // enf for tpi
	  //
	  Pop[p].Cell[i].Axonals=(Synapse *) calloc(Pop[p].Cell[i].Naxonals,sizeof(Synapse));

	  for(j=0;j<Pop[p].Cell[i].Naxonals;j++)
	    {
	      Pop[p].Cell[i].Axonals[j].TargetPop=tpni[j];
	      Pop[p].Cell[i].Axonals[j].TargetNeuron=ni[j];
              Pop[p].Cell[i].Axonals[j].TargetAxon=tpax[j];
              // check if the targeted axon exists
              if(tpax[j]!=-1){

                // check efficacy value
                if(eff[j]>1 || eff[j]<0){printf("For presynaptic inhibition, the efficacy should be between 0 and 1. Stop! \n"); exit(0);}

                flag_a=0;
                for(k=0;k<Pop[tpni[j]].Cell[i].Naxonals;k++){
                   if(PopD[tpni[j]].TargetPops[k]==tpax[j]){  // need to use PopD becuase the neuron may have not yet created.
                        Pop[p].Cell[i].Axonals[j].TargetAxon=k;
                        flag_a=1;
                     }
                 if(flag_a==1) break;
                  }

              }
		  Pop[p].Cell[i].Axonals[j].Ca_pre=1.0;
	      Pop[p].Cell[i].Axonals[j].Efficacy=eff[j];
	      Pop[p].Cell[i].Axonals[j].TargetReceptor=trni[j];
	      Pop[p].Cell[i].Axonals[j].D=D[j];
	      Pop[p].Cell[i].Axonals[j].pv=pv[j];
              Pop[p].Cell[i].Axonals[j].tauD=tauD[j];
	      Pop[p].Cell[i].Axonals[j].F=F[j];
	      Pop[p].Cell[i].Axonals[j].Fp=Fp[j];
              Pop[p].Cell[i].Axonals[j].tauF=tauF[j];
	      Pop[p].Cell[i].Axonals[j].LastConductance=ts[j];

              for(k=0;k<PopD[p].Nreceptors;k++){
                Pop[p].Cell[i].Axonals[j].St_pre[k]=0;
                Pop[p].Cell[i].Axonals[j].last_pre_update[k]=0;
                Pop[p].Cell[i].Axonals[j].Ca_pre_tau[k]=Pop[p].Cell[i].Tau[k];
              }
	    }
	} // end for i
    } // end for p


  if(flagverbose) { report("Network generated\n"); fflush(stdout); }
}


// same as generate network, but it does not allocate memory. Used in the multitrial mode
// (to change the network from trial to trial, the only way is to start again with a different seed)


int InitializeNetwork()
{

  int p,i,j,r,tp,tpi,tn,k,rd,na;
  int ni[MAXTN],tpni[MAXTN],trni[MAXTN];
  float eff[MAXTN],ts[MAXTN],timets[MAXTN];

  if(flagverbose) {
    report("\nGenerating network...\n");
  }

  dt=.1;
  delayindt=5;

  // loop over all the populations
  for(p=0;p<Npop;p++)
    {
      strcpy(Pop[p].Label,PopD[p].Label);
      if(flagverbose) {
	report("Population %2d: %s\n",p,Pop[p].Label);
      }

      for(rd=0;rd<PopD[p].Nreceptors;rd++)
	{
	  // external input: compute first the asymptoic mu and sigma (m,s in nphys notation) of the conductances
	  // CAREFUL: PROTOCOL FILE SHOULD BE REREAD TO INITIALIZE CORRECTLY EXT FREQS (THEY MIGHT HAVE BEEN CHANGED IN THE PREVIOUS TRIAL

	  PopD[p].MeanExtMuS[rd]=PopD[p].FreqExt[rd]*.001*PopD[p].MeanExtEff[rd]*PopD[p].MeanExtCon[rd]*PopD[p].Tau[rd];
	  PopD[p].MeanExtSigmaS[rd]=sqrt(PopD[p].Tau[rd]*.5*PopD[p].FreqExt[rd]*.001*PopD[p].MeanExtEff[rd]*PopD[p].MeanExtEff[rd]*PopD[p].MeanExtCon[rd]);
	  report("Pop %d Conductance: %d m=%f nS s=%f nS\n",p,rd,PopD[p].MeanExtMuS[rd],PopD[p].MeanExtSigmaS[rd]);
	}

      // init auxiliary variables like the tables of spikes
      if(delayindt>MAXDELAYINDT) { printf("ERROR: delay too long. Increase MAXDELAYINDT and recompile\n"); return -1;}
      Pop[p].CTableofSpikes=delayindt;
      Pop[p].DTableofSpikes=0;
      for(k=0;k<MAXDELAYINDT;k++)
	{
	  Pop[p].NTableofSpikes[k]=0;
	}

      // init     neurons and their axonal trees
      for(i=0;i<Pop[p].Ncells;i++)
	{
	  for(r=0;r<Pop[p].Cell[i].Nreceptors;r++)
	    {
	      Pop[p].Cell[i].LS[r]=0.;

	      // external input
	      Pop[p].Cell[i].ExtMuS[r]=PopD[p].MeanExtMuS[r];
	      Pop[p].Cell[i].ExtSigmaS[r]=PopD[p].MeanExtSigmaS[r];
	      Pop[p].Cell[i].ExtS[r]=0.;
	    } // end for r

	  // init
	  Pop[p].Cell[i].V=Pop[p].Cell[i].RestPot;
          Pop[p].Cell[i].Ca=0;
	  Pop[p].Cell[i].RefrState=0;
	  Pop[p].Cell[i].TimeLastSpike=-10000.; // time of the last emitted spike (for NMDA saturation)
	  Pop[p].Cell[i].PTimeLastSpike=-10000.; // time of spike emitted previous the last emitted spiek (for NMDA saturation)

	} // end for i
    } // end for p

  if(flagverbose) { report("Network initialized\n"); fflush(stdout); }
}



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

/*

Trick for NMDA:

$$ {ds_k \over dt} = -{s_k \over \tau} +\alpha(1-s_k)\sum_j \delta(t-t^k_j)$$

Multiply by $w_k$ and sum:

$$ {dS \over dt} = -S \over \tau+\sum_k \alpha (1-s_k)w_k \delta(t-t^k_j)$$

*/

// float current_freq;

int SimulateOneTimeStep()
{
  int aux;
  int p,i,j,r,sourceneuron,k;
  int tn,tp,tr,ta,tpi, tni, tri;
  float s,saturationfactor,total_St_pre;
  float Vaux; // auxiliary V: during the emission of the spike V is set artificially to 0. This is bad for the reversal potential
  float D,F,Ca_pre,DT;
  float g_adr, g_k, tau_max, alpha, beta, dv, n, tau_n, n_inif, efficacy, ExtMuS, ExtSigmaS;
  static float freq[MAXP][MAXRECEPTORS], mean_freq[MAXP][MAXRECEPTORS];
  static float last_freq[MAXP][MAXRECEPTORS];
  static int flag=0;
  static

  // DEBUG
  float nvalues=0,value=0;

  // END DEBUG

  if(flag==0){
    for(j=0;j<MAXP;j++){
    for(i=0;i<MAXRECEPTORS;i++)
      freq[j][i]=PopD[j].FreqExt[i];
      mean_freq[j][i]=PopD[j].FreqExt[i];
      last_freq[j][i]=PopD[j].FreqExt[i];
    }
   flag=1;
  }


  // Compute the decay of the total conductances and add external input
  // ------------------------------------------------------------------
  for(p=0;p<Npop;p++)
    {

      for(r=0;r<MAXRECEPTORS;r++){
	if(PopD[p].FreqExtSD[r]!=0){
	    // do{freq[p][r]=PopD[p].FreqExt[r]+gasdev()*PopD[p].FreqExtSD[r];}while(freq[p][r]<0);
	  do{freq[p][r]+=-(1-PopD[p].FreqExtDecay[r])*(last_freq[p][r]-PopD[p].FreqExt[r])
             +gasdev()*PopD[p].FreqExtSD[r];}while(freq[p][r]<0);
           last_freq[p][r]=freq[p][r]; ext_freq=freq[1][0];

	}
        else{
           if(PopD[p].FreqExtTau[r]==dt)
             {
//printf("cp1 PopD[%i].FreqExt[%i]=%f, FreqExtTau=%f \n",p,r,PopD[p].FreqExt[r],PopD[p].FreqExtTau[r]);
mean_freq[p][r]=PopD[p].FreqExt[r];}
           else{
//printf("cp2 PopD[%i].FreqExt[%i]=%f, FreqExtTau=%f, FreqExtDecay=%f \n",p,r,PopD[p].FreqExt[r],PopD[p].FreqExtTau[r],PopD[p].FreqExtDecay[r]);
             mean_freq[p][r]=(mean_freq[p][r]-PopD[p].FreqExt[r])*(1-dt/PopD[p].FreqExtTau[r])+PopD[p].FreqExt[r];}
             freq[p][r]=mean_freq[p][r]; ext_freq=freq[1][0];}
      }
      //      current_freq=freq[0][0];
      for(i=0;i<Pop[p].Ncells;i++)
	{
	  for(r=0;r<Pop[p].Cell[i].Nreceptors;r++)
	    {
              efficacy=PopD[p].MeanExtEff[r];
              Pop[p].Cell[i].ExtMuS[r]=freq[p][r]*.001*efficacy*PopD[p].MeanExtCon[r]*Pop[p].Cell[i].Tau[r];
              Pop[p].Cell[i].ExtSigmaS[r]=sqrt(Pop[p].Cell[i].Tau[r]*.5*freq[p][r]*.001*efficacy*efficacy*PopD[p].MeanExtCon[r]);
   	      s=Pop[p].Cell[i].ExtSigmaS[r];
	      if(s!=0.) // to optimize
		{
		  Pop[p].Cell[i].ExtS[r]+=dt/Pop[p].Cell[i].Tau[r]*(-Pop[p].Cell[i].ExtS[r]+Pop[p].Cell[i].ExtMuS[r])
		+s*sqrt(dt*2./Pop[p].Cell[i].Tau[r])*gauss();
		}
	      else {
		Pop[p].Cell[i].ExtS[r]+=dt/Pop[p].Cell[i].Tau[r]*(-Pop[p].Cell[i].ExtS[r]+Pop[p].Cell[i].ExtMuS[r]);
	      }

	      Pop[p].Cell[i].LS[r]*=exp(-dt/Pop[p].Cell[i].Tau[r]); // decay
	    }
	}
    }

  // Update the total conductances (changes provoked by the spikes)
  // --------------------------------------------------------------
  for(p=0;p<Npop;p++)
    {
      // loop over all the spikes emitted at time t-delay (they are received now)
      for(i=0;i<Pop[p].NTableofSpikes[Pop[p].DTableofSpikes];i++)
	{
	  sourceneuron=Pop[p].TableofSpikes[Pop[p].DTableofSpikes][i];

	  // for each spike, loop over the target conductances
	  for(j=0;j<Pop[p].Cell[sourceneuron].Naxonals;j++)
	    {
            if((ta=Pop[p].Cell[sourceneuron].Axonals[j].TargetAxon)!=-1){ //if this is a presynaptic inhibition
                tpi=Pop[p].Cell[sourceneuron].Axonals[j].TargetPop;
                tni=Pop[p].Cell[sourceneuron].Axonals[j].TargetNeuron;
                tri=Pop[p].Cell[sourceneuron].Axonals[j].TargetReceptor;
                if(DT=(Time-Pop[tpi].Cell[tni].Axonals[ta].last_pre_update[tri])>0){  // need to update presynaptic gating variable first
                   Pop[tpi].Cell[tni].Axonals[ta].St_pre[tri]=Pop[tpi].Cell[tni].Axonals[ta].St_pre[tri]*exp(-DT/Pop[tpi].Cell[tni].Axonals[ta].Ca_pre_tau[tri]);
                   Pop[tpi].Cell[tni].Axonals[ta].last_pre_update[tri]=Time;
                   }
                 // Now add the effect of presynaptic spikes
                 Pop[tpi].Cell[tni].Axonals[ta].St_pre[tri]+=Pop[p].Cell[sourceneuron].Axonals[j].Efficacy;
                 if(Pop[tpi].Cell[tni].Axonals[ta].St_pre[tri]>1) Pop[tpi].Cell[tni].Axonals[ta].St_pre[tri]=1.0;

              }
           else{ // Not a presynaptic inhibition, also need to update S_pre
               for(k=0;k<PopD[p].Nreceptors;k++){
                 if((strcmp(PopD[p].ReceptorLabel[k],"GABAA")==0)||(strcmp(PopD[p].ReceptorLabel[k],"GABAB")==0)){ // Only update GABAA and GABAB
                 DT=(Time-Pop[p].Cell[sourceneuron].Axonals[j].last_pre_update[k]);
                 if(DT>0&&(Pop[p].Cell[sourceneuron].Axonals[j].St_pre[k]>0)){
                    Pop[p].Cell[sourceneuron].Axonals[j].St_pre[k]=Pop[p].Cell[sourceneuron].Axonals[j].St_pre[k]*exp(-DT/Pop[p].Cell[sourceneuron].Axonals[j].Ca_pre_tau[k]);
                    Pop[p].Cell[sourceneuron].Axonals[j].last_pre_update[k]=Time;
                      }
                  }
                 }

             // Short-term depression: Calculate the recovery of D first.

               Pop[p].Cell[sourceneuron].Axonals[j].D = 1-(1-Pop[p].Cell[sourceneuron].Axonals[j].D)*exp(-(Time-Pop[p].Cell[sourceneuron].PTimeLastSpike)/Pop[p].Cell[sourceneuron].Axonals[j].tauD);

	        D=Pop[p].Cell[sourceneuron].Axonals[j].D;

               total_St_pre=0;
               for(k=1;k<PopD[p].Nreceptors;k++){
                  total_St_pre+=Pop[p].Cell[sourceneuron].Axonals[j].St_pre[k];
                 }
              Ca_pre=1-total_St_pre;
             //D=1.0;
	     // if(sourceneuron==20) printf(" %.1f %.1f .1%f .3%f \n",Time, Pop[p].Cell[sourceneuron].PTimeLastSpike, Pop[p].Cell[sourceneuron].tauD, Pop[p].Cell[sourceneuron].D);

             // Short-term facilitation: Calculate the F value
              if(Pop[p].Cell[sourceneuron].Axonals[j].Fp==0)
                F=1;
              else{
               Pop[p].Cell[sourceneuron].Axonals[j].F = Pop[p].Cell[sourceneuron].Axonals[j].F*exp(-(Time-Pop[p].Cell[sourceneuron].PTimeLastSpike)/Pop[p].Cell[sourceneuron].Axonals[j].tauF);
               F=Pop[p].Cell[sourceneuron].Axonals[j].F;
	        }
	       tn=Pop[p].Cell[sourceneuron].Axonals[j].TargetNeuron;
	       tp=Pop[p].Cell[sourceneuron].Axonals[j].TargetPop;
	       tr=Pop[p].Cell[sourceneuron].Axonals[j].TargetReceptor;

	       if(Pop[p].Cell[sourceneuron].Axonals[j].LastConductance<0.) { // NO NMDA (no saturation)
		  Pop[tp].Cell[tn].LS[tr]+=powf(Ca_pre,3.5)*D*F*Pop[p].Cell[sourceneuron].Axonals[j].Efficacy; // no saturation
	        }
	       else {

		// Now it should be correct. ALPHA factor to be determined (jump for every spike): now it is the maximum of a single spike
		// TEMPORARY
#define ALPHA 0.6332
		// 0.6332 is the best value for the area at 55 Hz. Best value depends on the frequency!!!
		 Pop[p].Cell[sourceneuron].Axonals[j].LastConductance*=exp(-(Time-Pop[p].Cell[sourceneuron].PTimeLastSpike)/Pop[tp].Cell[tn].Tau[tr]);

		 Pop[tp].Cell[tn].LS[tr]+=D*F*Pop[p].Cell[sourceneuron].Axonals[j].Efficacy*(1.-Pop[p].Cell[sourceneuron].Axonals[j].LastConductance)*ALPHA;

		 Pop[p].Cell[sourceneuron].Axonals[j].LastConductance+=ALPHA*(1.-Pop[p].Cell[sourceneuron].Axonals[j].LastConductance);

	       }
 // Calculate the change of D and F due to the spike.
                Pop[p].Cell[sourceneuron].Axonals[j].F += Pop[p].Cell[sourceneuron].Axonals[j].Fp*(1-Pop[p].Cell[sourceneuron].Axonals[j].F);

               Pop[p].Cell[sourceneuron].Axonals[j].D -= Pop[p].Cell[sourceneuron].Axonals[j].pv*powf(Ca_pre,3.5)*F*Pop[p].Cell[sourceneuron].Axonals[j].D;
                if(Pop[p].Cell[sourceneuron].Axonals[j].D<0)  Pop[p].Cell[sourceneuron].Axonals[j].D=0;


              }
	    } // for j
		//                if(sourceneuron==20) printf("   %f\n",Pop[p].Cell[sourceneuron].D);
                //if(Pop[p].Cell[sourceneuron].D<=0) Pop[p].Cell[sourceneuron].D=0;
	}
    }

  // Update the neuronal variables
  // -----------------------------

  for(p=0;p<Npop;p++)
    {
      Pop[p].NTableofSpikes[Pop[p].CTableofSpikes]=0; // reset the number of emitted spikes for pop p
      for(i=0;i<Pop[p].Ncells;i++)
	{

	  if(Pop[p].Cell[i].V>Pop[p].Cell[i].Threshold)
	    {
	      Pop[p].Cell[i].V=Pop[p].Cell[i].ResetPot; // special state after emission
	      Pop[p].Cell[i].RefrState--;
	      continue;
	    }

	  // refractory period

	  if(Pop[p].Cell[i].RefrState) {
	    Pop[p].Cell[i].RefrState--;
	    continue;
	  }

          //Anomalous delayed rectiflier (ADR)
          if(Pop[p].Cell[i].g_adr_max==0) // No ADR
             g_adr=0;
          else //with ADR
	     g_adr = Pop[p].Cell[i].g_adr_max/(1+exp((Pop[p].Cell[i].V-Pop[p].Cell[i].Vadr_h)/Pop[p].Cell[i].Vadr_s));

          //potassium outward rectifying current
          if(Pop[p].Cell[i].g_k_max==0) // No outward current
             g_k=0;
          else{ // with outward current
	    //  g_k = Pop[p].Cell[i].g_k_max/(1+exp((-Pop[p].Cell[i].V+Pop[p].Cell[i].Vk_h)/Pop[p].Cell[i].Vk_s));
            tau_max=Pop[p].Cell[i].tau_k_max;
            dv=(Pop[p].Cell[i].V+55.0);
	    //  alpha=(-10.0/tau_max)*(dv-49)/(1-exp(-(dv-49)/100));
	    //            beta=(0.17/tau_max)*exp(-(dv-11)/10);
            tau_n=tau_max/(exp(-1*dv/30)+exp(dv/30));
            n_inif=1/(1+exp(-(Pop[p].Cell[i].V-Pop[p].Cell[i].Vk_h)/Pop[p].Cell[i].Vk_s));

            Pop[p].Cell[i].n_k+=-1/tau_n*(Pop[p].Cell[i].n_k-n_inif);
            g_k=Pop[p].Cell[i].g_k_max*Pop[p].Cell[i].n_k;
	    //	    printf("v=%f dv=%f alpha=%f beta=%f g_k=%f \n",Pop[p].Cell[i].V,dv,alpha,beta,g_k);
	  }


          // decay
          if(Pop[p].Cell[i].g_ahp==0) //No spike-frequency adaptation
    	    Pop[p].Cell[i].V+= -dt*(1/Pop[p].Cell[i].Taum*(Pop[p].Cell[i].V-Pop[p].Cell[i].RestPot)
                               + g_adr/Pop[p].Cell[i].C*(Pop[p].Cell[i].V-Pop[p].Cell[i].ADRRevPot)
                               + g_k/Pop[p].Cell[i].C*(Pop[p].Cell[i].V-Pop[p].Cell[i].ADRRevPot));
          else // with spike-frequency adaptation (the factor 1/1000 is needed to convert nS/nF to 1/ms)
     	    Pop[p].Cell[i].V+= -dt*(1/Pop[p].Cell[i].Taum *(Pop[p].Cell[i].V-Pop[p].Cell[i].RestPot)
                + Pop[p].Cell[i].Ca*Pop[p].Cell[i].g_ahp/Pop[p].Cell[i].C*0.001*(Pop[p].Cell[i].V-Pop[p].Cell[i].Vk)
                + g_adr/Pop[p].Cell[i].C*(Pop[p].Cell[i].V-Pop[p].Cell[i].ADRRevPot)
                + g_k/Pop[p].Cell[i].C*(Pop[p].Cell[i].V-Pop[p].Cell[i].ADRRevPot));


          // [Ca] decay -- 1st order approximation
          Pop[p].Cell[i].Ca += -Pop[p].Cell[i].Ca*dt/Pop[p].Cell[i].tau_ca;

	  Vaux=Pop[p].Cell[i].V;
      	  if(Vaux>Pop[p].Cell[i].Threshold) Vaux=Pop[p].Cell[i].Threshold;

	  // now add the synaptic currents (the factor 1/1000 is needed to convert nS/nF to 1/ms)
	  for(r=0;r<Pop[p].Cell[i].Nreceptors;r++)
	    {
	      if(Pop[p].Cell[i].MgFlag[r]) { // magnesium block
		Pop[p].Cell[i].V+=dt*(Pop[p].Cell[i].RevPots[r]-Vaux)*.001*(Pop[p].Cell[i].LS[r]+Pop[p].Cell[i].ExtS[r])
		  /Pop[p].Cell[i].C/(1.+exp(-0.062*Vaux)/3.57);

		/*		// DEBUG DEBUG DEBUG
		if(p==1 && i==0) {
		  printf("[%f]",Pop[p].Cell[i].LS[r]);
		}
		// end DEBUG */

	      } else {
	      Pop[p].Cell[i].V+=dt*(Pop[p].Cell[i].RevPots[r]-Vaux)*.001*(Pop[p].Cell[i].LS[r]+Pop[p].Cell[i].ExtS[r])
              /Pop[p].Cell[i].C;
	      }
	    }
	  // spike condition
	  if(Pop[p].Cell[i].V>Pop[p].Cell[i].Threshold) {

	    // a spike is emitted
	    Pop[p].TableofSpikes[Pop[p].CTableofSpikes][Pop[p].NTableofSpikes[Pop[p].CTableofSpikes]]=i;
	    if(Pop[p].NTableofSpikes[Pop[p].CTableofSpikes]<MAXSPIKESINDT-1) Pop[p].NTableofSpikes[Pop[p].CTableofSpikes]++;
	    else { printf("\nERROR: too many spikes in a dt (change MAXSPIKESINDT and recompile)\n"); fflush(stdout); }

	    Pop[p].Cell[i].V=Pop[p].Cell[i].ResetPot;
	    Pop[p].Cell[i].PTimeLastSpike=Pop[p].Cell[i].TimeLastSpike;
	    Pop[p].Cell[i].TimeLastSpike=Time;

	    Pop[p].Cell[i].V=0.; // spike! (temporary)
	    Pop[p].Cell[i].RefrState=(int)floor(Pop[p].Cell[i].RefractT/dt); // refractory period!!! (temporary)

            // [Ca] increases;
            Pop[p].Cell[i].Ca += Pop[p].Cell[i].alpha_ca;

	  }
	}
    }

  // Update the pointers of the table of spikes
  // ---------------------------------------------------------------

  for(p=0;p<Npop;p++)
    {
      Pop[p].CTableofSpikes++; if(Pop[p].CTableofSpikes>MAXDELAYINDT-1) Pop[p].CTableofSpikes=0;
      Pop[p].DTableofSpikes++; if(Pop[p].DTableofSpikes>MAXDELAYINDT-1) Pop[p].DTableofSpikes=0;
    }
}

// DATA ANALYSIS
// --------------------------------------------------------------------------------------------

// A matlab script is generated at the beginning, just before starting the simulation.

#define TIMEWINDOWFORFREQ 20.  // time window on which the mean pop frequency is estimated (in ms)
#define RASTERPLOTNEURONS 50   // number of neurons in the raster plot
#define NUMBEROFTRACES 2  // number of visualized traces of V
#define STEPSFORPRINTIGFREQS 500 // every ... step, the frequencies are printed
#define STEPSFORSAVINGFREQS 10 // every ... step the frequencies are saved (not always, to save space)
#define STEPSFORFLUSHING 80
#define SPIKEBUFFER 300 //szie of the array SpikeBuffer[] in SaveSpikes(). This value should be no
                        // smaller than TIMEWINDOWFORFREQ/dt

int NumberofTraces=5;


int GenMatLabFile()
{
  int i,j,nsp,r;
  FILE *devmatlab;
  static char color[5]="rgby";

  devmatlab=fopen("mrealsimu.m","w");
  if(devmatlab==NULL) { printf("ERROR: Unable to open mathlab file\n"); return 0;}

  fprintf(devmatlab,"clear all\n");
  fprintf(devmatlab,"load popfreqs.dat\n");
  if(NumberofTraces>0) {
    fprintf(devmatlab,"load poptraces.dat\n");
  }

  nsp=2+2*NumberofTraces; // total number of subplots

  for(i=0;i<Npop;i++)
    {
      // plot the raster
      fprintf(devmatlab,"figure(%d)\nclf\nsubplot(%d,1,1);\n",i+1,nsp);
      fprintf(devmatlab,"load pop%d.dat\n",i+1);
      fprintf(devmatlab,"plotRast(pop%d,0.,%f,%f,%d);\n",i+1,TrialDuration,dt,RASTERPLOTNEURONS);
      fprintf(devmatlab,"title('%s - Raster plot');\n",Pop[i].Label);
      fprintf(devmatlab,"xlabel('Time (ms)');\n");

      // plot the mean frequencies
      fprintf(devmatlab,"subplot(%d,1,2);\nplot(popfreqs(:,1),popfreqs(:,%d));\n",nsp,i+2);
      fprintf(devmatlab,"ylabel('Mean spike freq (Hz)');\nxlabel('Time (ms)');\n");
      fprintf(devmatlab,"set(gca,'XLim',[0 %f]);\n",TrialDuration);

      // plot the traces
      for(j=0;j<NumberofTraces;j++)
	{
	  fprintf(devmatlab,"subplot(%d,1,%d);\n",nsp,3+j*2);
	  fprintf(devmatlab,"plot(poptraces(:,1),poptraces(:,%d));\n",2+j*(Pop[i].Cell[j].Nreceptors+1)+i*(NumberofTraces*(Pop[i].Cell[j].Nreceptors+1)));

	  fprintf(devmatlab,"ylabel('Membrane pot (mV) - Neuron %d')\n",j);
	  fprintf(devmatlab,"set(gca,'YLim',[%f 5]);\nset(gca,'XLim',[0 %f]);\n",PopD[i].RestPot,TrialDuration);
	  fprintf(devmatlab,"subplot(%d,1,%d);\n",nsp,3+j*2+1);

	  for(r=0;r<Pop[i].Cell[j].Nreceptors;r++)
	    {
	      fprintf(devmatlab,"plot(poptraces(:,1),poptraces(:,%d),'%c');\nhold on;\n",r+3+j*(Pop[i].Cell[j].Nreceptors+1)+i*(NumberofTraces*(Pop[i].Cell[j].Nreceptors+1)),color[r]);
	    }
	  fprintf(devmatlab,"set(gca,'XLim',[0 %f]);\n",TrialDuration);
	  fprintf(devmatlab,"title('Tot conductances (nS) - Neuron %d')\n",j);
	}
    }

  fclose(devmatlab);

  return 1;
}


int GenMatLabMultiTrial()
{
  FILE *devmatlab;
  int p,trial;

  devmatlab=fopen("mrealmulti.m","w");
  if(devmatlab==NULL) { printf("ERROR: Unable to open multi trial mathlab file\n"); return 0;}

  fprintf(devmatlab,"clear all\nfigure(1);\nclf\n");

  // load all frequency files
  for(trial=0;trial<CurrentTrial+1;trial++)
    {
      fprintf(devmatlab,"load popfreqs%d.dat;\n",trial);
    }

  for(p=0;p<Npop;p++)
    {
      fprintf(devmatlab,"subplot(%d,1,%d);\n",Npop,p+1);
      for(trial=0;trial<CurrentTrial+1;trial++)
	{
	     fprintf(devmatlab,"plot(popfreqs%d(:,1),popfreqs%d(:,%d));\nhold on\n",trial,trial,p+2);
	}
      fprintf(devmatlab,"ylabel('Mean spike freq (Hz)');\nxlabel('Time (ms)');\n");
      fprintf(devmatlab,"set(gca,'XLim',[0 %f]);\n",TrialDuration);
      fprintf(devmatlab,"title('Population: %s');\n",Pop[p].Label);
    }
  fclose(devmatlab);
}

int SaveSpikes(int eventflag, long sno)
{

  static int InitFlag=1;
  static FILE *devspikes[MAXP],*devfreqs;
  static float meanfreqs[MAXP];
  static float timelastevent;
  static float meanfreqsbetweenevents[MAXP];
  static float SpikeBuffer[MAXP][SPIKEBUFFER];
  static int counter;
  static int lasttrial=0;
  static int currentpt, buffersize;  // for SpikeBuffer[][]
  int i,p,TempPointer;
  float  TempBuffer;
  char TempName[100];

  // initialize if it is the first call
  if(InitFlag || (lasttrial!=CurrentTrial)) {

    // if there is a new trial, close all the files of the previous trial
    if(lasttrial!=CurrentTrial) {
      if(FlagSaveAllSpikes)
	{
	  for(p=0;p<Npop;p++)
	    {
	      fclose(devspikes[p]);
	    }
	}
      lasttrial=CurrentTrial;
      // and then open th new ones

    }

    buffersize=(int)(TIMEWINDOWFORFREQ/dt);
    if(buffersize>=SPIKEBUFFER) {printf("Error: SPIKEBUFFER is too small. Edit the code and recomplie the program. \n"); exit(1);};
    currentpt=0;

    // open all files
    printf("\n Time (ms) ");

    for(p=0;p<Npop;p++)
      {
	sprintf(TempName,"pop%d_%d.%ld.dat",p+1,CurrentTrial,sno);
	if(FlagSaveAllSpikes)
	  {
	    devspikes[p]=fopen(TempName,"w");
	    if(devspikes[p]==NULL) return 0;
	  }
	meanfreqs[p]=0.;
	meanfreqsbetweenevents[p]=0.;

	for(i=0;i<SPIKEBUFFER;i++)
	  SpikeBuffer[p][i]=0;

	printf("%12s",Pop[p].Label);
      }
    timelastevent=0.;

    sprintf(TempName,"popfreqs%d.%ld.dat",CurrentTrial,sno);
    devfreqs=fopen(TempName,"w");
    if(devfreqs==NULL) return 0;
    InitFlag=0;
    counter=0;
    printf("\n---------------------------------------------------------------\n");
  }

  if((counter % STEPSFORSAVINGFREQS)==0) fprintf(devfreqs,"%f\t",Time);

  if((counter % STEPSFORPRINTIGFREQS)==0) printf("%7.1f ms",Time);

  for(p=0;p<Npop;p++)
    {
      TempPointer=Pop[p].CTableofSpikes-1;
      if(TempPointer<0) TempPointer=MAXDELAYINDT-1;

      meanfreqsbetweenevents[p]+=(float)Pop[p].NTableofSpikes[TempPointer]/(float)Pop[p].Ncells/dt*1000.;


      meanfreqs[p] += -SpikeBuffer[p][currentpt] + (TempBuffer=(float)Pop[p].NTableofSpikes[TempPointer]/(float)Pop[p].Ncells*1000/TIMEWINDOWFORFREQ);
      SpikeBuffer[p][currentpt] = TempBuffer;

      if(currentpt<buffersize-1)currentpt++; else currentpt=0;
      //      meanfreqs[p]+=dt/TIMEWINDOWFORFREQ*(-meanfreqs[p]+(float)Pop[p].NTableofSpikes[TempPointer]/(float)Pop[p].Ncells/dt*1000.); // compute mean freq in Hz on a time window of 10 ms

      if((counter % STEPSFORPRINTIGFREQS)==0) printf("%12.1f",meanfreqs[p]);
      //      if((counter % STEPSFORPRINTIGFREQS)==0) printf("%s FreqExtSD=%f \n",PopD[0].ReceptorLabel[0],PopD[0].FreqExtSD[0]);
      if((counter % STEPSFORSAVINGFREQS)==0) fprintf(devfreqs,"%f\t",meanfreqs[p]); // mean frequency in Hz per neuron

      if(FlagSaveAllSpikes)
	{
	  for(i=0;i<Pop[p].NTableofSpikes[TempPointer];i++) {
	    fprintf(devspikes[p],"%d\t%f\n",Pop[p].TableofSpikes[TempPointer][i],Time);
	  }
	}
    }

    if((counter % STEPSFORSAVINGFREQS)==0) fprintf(devfreqs,"\n");
    if((counter % STEPSFORPRINTIGFREQS)==0) printf("\n");

  //  if((counter % STEPSFORSAVINGFREQS)==0) fprintf(devfreqs," %f \n",current_freq);
  //  if((counter % STEPSFORPRINTIGFREQS)==0) printf(" %f \n",current_freq);

  if((counter % STEPSFORFLUSHING)==0) {
    if(FlagSaveAllSpikes)
      {
	for(p=0;p<Npop;p++)
	  {
	    fflush(devspikes[p]);
	  }
      }
    fflush(devfreqs);
  }

  if(eventflag) {

    printf("Average:  ");
    for(p=0;p<Npop;p++)
      {
	printf("%12.1f",meanfreqsbetweenevents[p]/(Time-timelastevent)*dt);
	meanfreqsbetweenevents[p]=0.;
      }
    printf("\n");
    fflush(stdout);
    timelastevent=Time;
  }

  counter++;

  return 1;
}

int SaveTraces(long sno)
{
  int i,p,r;
  static int FlagInit=1;
  static FILE *devtraces;
  static int lasttrial=0;
  char TempName[100];
  float Vaux;

  if(NumberofTraces==0) return 1;

  if(FlagInit || (lasttrial!=CurrentTrial)) {
    if(lasttrial!=CurrentTrial) {
      fclose(devtraces);
      lasttrial=CurrentTrial;
    }

    sprintf(TempName,"poptraces%d.%ld.dat",CurrentTrial,sno);
    devtraces=fopen(TempName,"w");
    if(devtraces==NULL) return 0;
    FlagInit=0;
  }

  fprintf(devtraces,"%f\t",Time);
  for(p=0;p<Npop;p++)
    {
      for(i=0;i<NumberofTraces;i++)
	{
	  Vaux=Pop[p].Cell[i].V;
	  if(Vaux>Pop[p].Cell[i].Threshold) Vaux=Pop[p].Cell[i].Threshold;

	  fprintf(devtraces,"%f\t",Pop[p].Cell[i].V);
	  for(r=0;r<Pop[p].Cell[i].Nreceptors;r++)
	    {
	      if(Pop[p].Cell[i].MgFlag[r]) {
		fprintf(devtraces,"%f\t",Pop[p].Cell[i].LS[r]);  // /(1.+exp(-0.062*Vaux)/3.57)); mg block now not included
	      }
	      else {
		fprintf(devtraces,"%f\t",Pop[p].Cell[i].LS[r]);
	      }
	    }
	}
    }
  fprintf(devtraces,"\n");
}


// Handle event
// -------------------------------------------------------------------------------------------------------------------

int HandleEvent(void)
{
  int i,p,r,q,j;
  float efficacy,MeanEff;

  if(Events[CEvent].Type==ENDOFTRIAL) {
    return 0;
  }

  if(Events[CEvent].Type==CHANGEEXTFREQ) {
    p=Events[CEvent].PopNumber;
    r=Events[CEvent].ReceptorNumber;
    PopD[p].FreqExt[r]=Events[CEvent].FreqExt;
    if(Events[CEvent].FreqExtsd!=0)
       PopD[p].FreqExt[r]=Events[CEvent].FreqExt+gasdev()*Events[CEvent].FreqExtsd;
    else
       PopD[p].FreqExt[r]=Events[CEvent].FreqExt;

 //   printf("Events[CEvent].FreqExtDecay=%f \n",Events[CEvent].FreqExtDecay);
    PopD[p].FreqExtTau[r]=Events[CEvent].FreqExtDecay;
    printf("%7.1f ------------------ Event: %s ----------------\n",Time,Events[CEvent].Label);

    efficacy=PopD[p].MeanExtEff[r];
    PopD[p].MeanExtMuS[r]=PopD[p].FreqExt[r]*.001*efficacy*PopD[p].MeanExtCon[r]*PopD[p].Tau[r];
    PopD[p].MeanExtSigmaS[r]=sqrt(PopD[p].Tau[r]*.5*PopD[p].FreqExt[r]*.001*efficacy*efficacy*PopD[p].MeanExtCon[r]);

    for(i=0;i<Pop[p].Ncells;i++)
      {
        do{efficacy=(gasdev()*PopD[p].ExtEffSD[r])+PopD[p].MeanExtEff[r];}while(efficacy<0);
	Pop[p].Cell[i].ExtMuS[r]=PopD[p].FreqExt[r]*.001*efficacy*PopD[p].MeanExtCon[r]*PopD[p].Tau[r];
	Pop[p].Cell[i].ExtSigmaS[r]=sqrt(PopD[p].Tau[r]*.5*PopD[p].FreqExt[r]*.001*efficacy*efficacy*PopD[p].MeanExtCon[r]);
      }
  }

  if(Events[CEvent].Type==CHANGEEXTFREQSD) {
    p=Events[CEvent].PopNumber;
    r=Events[CEvent].ReceptorNumber;
    PopD[p].FreqExtSD[r]=Events[CEvent].FreqExtSD;
    printf("%7.1f ------------------ Event: %s ----------------\n",Time,Events[CEvent].Label);


  }

 if(Events[CEvent].Type==CHANGEMEANEFF) {
    p=Events[CEvent].PopNumber;
    q=Events[CEvent].TargetPopNumber;
    r=Events[CEvent].ReceptorNumber;
    printf("%7.1f ------------------ Event: %s ----------------\n",Time,Events[CEvent].Label);
    MeanEff=Events[CEvent].MeanEff;
    for(i=0;i<PopD[p].NTargetPops;i++)
       if(PopD[p].TargetPops[i]==q){
	 if(PopD[p].SynP[i].TargetReceptor=r)
           PopD[p].SynP[i].MeanEfficacy=MeanEff;
       }
    for(i=0;i<Pop[p].Ncells;i++){
      for(j=0;j<Pop[p].Cell[i].Naxonals;j++)
	if((Pop[p].Cell[i].Axonals[j].TargetPop==q)&&(Pop[p].Cell[i].Axonals[j].TargetReceptor==r))
          Pop[p].Cell[i].Axonals[j].Efficacy=MeanEff;
    }
 }

  CEvent++;
  NextEventTime=Events[CEvent].ETime;

  return 1;
}


main(int argc, char *argv[])
{
  int ti,runflag,tistepforsaving,rseed;
  long iseed;
  long sno=0;

  iseed=1;
  sran1(&iseed);

  //  char prefix[50];

  // read line commands

  flagverbose=0;
  NumberofTrials=1;
  FlagSaveAllSpikes=1;
  //  strcpy(prefix,"cl_");

  if(argc>1) {
    do {
      if(strncmp(argv[argc-1],"-v",2)==0) { flagverbose=1; argc--; continue; }
      if(strncmp(argv[argc-1],"-h",2)==0) {
	printf("realsimu - Ver. 0.8\n Usage:\n-h  : this help\n-v  : verbose mode\n-t# : number of saved traces per population\n");
	printf("-T# : number of trials (the network is the same for each trial, the realization of the ext noise changes)\n");
	printf("-ns : spikes and traces are not saved. Only the mean frequencies are saved for each trial\n");
        printf("-s# : seed for the random number generator.\n");
	return -1;
      }
      if(strncmp(argv[argc-1],"-t",2)==0) {
	NumberofTraces=atoi(&argv[argc-1][2]);
	printf("Number of saved traces: %d\n",NumberofTraces);
	argc--; continue;
      }
      if(strncmp(argv[argc-1],"-s",2)==0) {
	rseed=atoi(&argv[argc-1][2]);
	srand49(rseed);
        iseed=-1*(long)rseed;
        sran1(&iseed);
	printf("Seed for random generator: %d\n",rseed);
	argc--; continue;
      }
      if(strncmp(argv[argc-1],"-ns",3)==0) {
	FlagSaveAllSpikes=0;
	NumberofTraces=0;
	printf("Spikes are not saved\n");
	argc--; continue;
      }
      if(strncmp(argv[argc-1],"-no",3)==0) {
	sno=atol(&argv[argc-1][3]);
	printf("Output file sequence: %ld\n",sno);
	argc--; continue;
      }
      if(strncmp(argv[argc-1],"-T",2)==0) {
	NumberofTrials=atoi(&argv[argc-1][2]);
	printf("Number of trials: %d\n",NumberofTrials);
	argc--; continue;
      }
      /*
      if(strncmp(argv[argc-1],"-p",2)==0) {
	strcpy=(prefix,&argv[argc-1][2]);
	printf("Prefix to all input/output files: %s\n",prefix);
	argc--; continue;
      }
      */
      printf("ERROR: unrecognized option: %s\n",argv[argc]);
      argc--;
    } while(argc>1);
  }

  if(DescribeNetwork()==-1) {printf("Unrecoverable error in parsing the conf file, exiting...\n"); return -1;}
  if(ParseProtocol()==-1) { printf("Unrecoverable error in parsing the protocol file, exiting...\n"); return -1;}

  GenerateNetwork();
  GenMatLabFile();

  for(CurrentTrial=0;CurrentTrial<NumberofTrials;CurrentTrial++)
    {
      report("Trial #\%d\n=====================================================================\n",CurrentTrial);
      CEvent=0;
      NextEventTime=Events[CEvent].ETime;
      runflag=1;

      GenMatLabMultiTrial();

      if(CurrentTrial) InitializeNetwork();

      for(ti=0,Time=0.;runflag;Time+=dt,ti++)
	{
	  SimulateOneTimeStep();
	  // handle all the events
	  if(Time>=NextEventTime) {
	    SaveSpikes(1,sno);
	    while(Time>=NextEventTime)
	      {
		runflag=HandleEvent();
		if(runflag==0) break;
	      }
	  }
	  else SaveSpikes(0,sno);

	  SaveTraces(sno);
	}
      report("\rEnd of the trial\n");
    }
}

Loading data, please wait...