Studies of stimulus parameters for seizure disruption using NN simulations (Anderson et al. 2007)

 Download zip file 
Help downloading and running models
Accession:98902
Architecturally realistic neocortical model using seven classes of excitatory and inhibitory single compartment Hodgkin-Huxley cells. Wiring is adapted to minicolumn hypothesis and incorporates visual and neocortical data. Simulation demonstrates spontaneous bursting onset and cessation, and activity can be altered with external electric field.
Reference:
1 . Anderson WS, Kudela P, Cho J, Bergey GK, Franaszczuk PJ (2007) Studies of stimulus parameters for seizure disruption using neural network simulations. Biol Cybern 97:173-94 [PubMed]
Citations  Citation Browser
Model Information (Click on a link to find other models with that property)
Model Type: Realistic Network;
Brain Region(s)/Organism:
Cell Type(s): Neocortex L5/6 pyramidal GLU cell; Neocortex L2/3 pyramidal GLU cell; Neocortex V1 interneuron basket PV GABA cell; Neocortex fast spiking (FS) interneuron;
Channel(s): I A; I K; I K,leak; I K,Ca; I Sodium; I Calcium;
Gap Junctions:
Receptor(s): AMPA; Gaba;
Gene(s):
Transmitter(s):
Simulation Environment: C or C++ program;
Model Concept(s): Bursting; Epilepsy; Vision;
Implementer(s):
Search NeuronDB for information about:  Neocortex L5/6 pyramidal GLU cell; Neocortex L2/3 pyramidal GLU cell; Neocortex V1 interneuron basket PV GABA cell; AMPA; Gaba; I A; I K; I K,leak; I K,Ca; I Sodium; I Calcium;
/* Program do generacji polaczen synaptycznych 
   na wyjsciu ma produkowac plik z polaczeniami z zapisana para neuronow waga i
   delay. Nastenamepny program (lub druga czesc tego samego, moze wykorzystac to
   do generacji struktury neuronow i synaps we wlasciwym formacie.

   05/29/2003 P. Kudela, now multi-class supported  
   07/27/2005 S. Anderson 
*/
#define RISC
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sys/types.h>
#include <time.h>
#include <math.h>
#include "lnet.h"

#define randomm(num) (int)((double)num*(double)random()/((double)(RAND_MAX+1.0)))
#define randomize()     srandom((unsigned)time(NULL))

int lrodz,histoint;               /* number of neuron classes */
int *NN;                 /* NN[lrodz] total number of neurons within each class */
int *NOF;                /* NOF[lrodz] table of index offsets */ 
int *NB; 

int **N;             /* N[lrodz][3] size of subnetwork for each class   */
int **k;             /* k[lrodz][3] scaling factor to align class networks  */

int ***o;            /* o[lrodz][lrodz][3] neighborhood array  */
int **s;             /* s[lrodz][lrodz] synapsy na we/wy z otoczenie */
int **w,**sdw;    /* w[lrodz][lrodz] sdw[][] synaptic weights and dispersions */
int **d,**sdd;    /* d[lrodz][lrodz] sdd[][] synaptic delays and dispersions */

int *tab;         /* temporary neighborhood array tab[maxot]*/ 

int get_kind(int nr)
{
  int kind=0;
  while(NOF[++kind]<=nr);
  return kind-1;
}

struct SYN_TMP
{
NEURON_INDEX    neuron_pre;
NEURON_INDEX    neuron_post;
NEURON_INDEX    preclass_num;
WEIGHT          weight;
DELAY           delay;
};

int sort_delay(const void *, const void *);
int sort_neu_num(const void *, const void *);
void set_neu(struct NEURON *);
unsigned seed_save();

FILE *cfg;
int histo_f,out_f,vout_f,noise_f,iext_f,link_f,stim_f;  /* flags for simulator */
char flags[100];

void  set_flags()
{
  /*signed*/ char z;
  /*unsigned*/ char buf[100],*l;
  l=flags;
  histo_f=out_f=vout_f=noise_f=iext_f=link_f=stim_f=0;  /* flags for simulator */

  fscanf(cfg,"%s",flags);
  do
    {
      sscanf(l,"%c,",&z);
      l+=2;

      switch(z)
	{
	case 'N': noise_f = 1;
	  break;
	case 'H': histo_f = 1;
	  break;
	case 'V': vout_f = 1;
	  break;
	case 'O': out_f = 1;
	  break;
	case 'I': iext_f = 1;
	  break;
	case 'L': link_f = 1;
	  break;
	case 'S': stim_f = 1;
	  iext_f = 1;
	case '$': break;
	default :
	  fprintf(stderr,"Unknown option %c in cfg",z);
	}
    }while (z != '$');

	  
  strcpy(buf,"Options:");
  if (noise_f) strcat(buf,"NOISE,");
  if (iext_f) strcat(buf,"IEXT,");
  if (histo_f) strcat(buf,"HISTO,");
  if (out_f) strcat(buf,"OUT,");
  if (vout_f) strcat(buf,"VOUT,");
  if (link_f) strcat(buf,"LINK,");
  if (stim_f) strcat(buf,"STIM,");
  puts(buf);
	  
}

int main(int argc, char  *argv[])
{
 int i,j,l,m,n,jp;
 int nr,ir,ix,iy,iz;  /* index for neurons */
 int maxot, tmpot;
 int xyz[3]={0,0,0};
 int no,nl;
 ulong ns,is;
 int *neu_pre;
 int Nneu;
 int torus_flag;
 int syn_n;
 char label[]="xyz",name[100],fname[100],buff[10],strdummy[10];
 FILE *inputfile,*outputfile,*par,*synfile,*neufile,*opar,*nodefile;
 struct SYN_TMP *syn_tmp;
 struct SYNAPSE syn;
 struct NEURON *neu;
 struct PARAMETERS p;
 int nproc,nsub,iproc,isub,vsub,nproc2,nsub2;
 SPIKE_INTERVAL delays[MAX_KIND][MAX_KIND];
 int otoczenie(int,int,int*,int);
 int otoczenie2(int,int,int*,int);
 int losuj(int*,int);
 int dummy,gumba;
 double testrandnoise;
 extern int zapisz(int,int,int,int,FILE *);
 extern int sort_neu_num(const void *, const void *);

 int **flag;
 int **flagv; 
 if(argc<5)
   {fprintf(stderr,"USAGE:\n\tmknodes <inp_name> <out_fname> <sqrt(subnets/node)> <sqrt(nodes)> [<seed>]\n");
   return 1;
   }

  if(argc>5)
   srandom((unsigned)atoi(argv[5]));
  else
   printf("seed %u\n",seed_save()); /* seed for srandom */

 strcpy(name,argv[1]);
  inputfile = fopen(strcat(name,".sub"),"r");
  if(inputfile==NULL){fprintf(stderr, "Subnetwork definition file not present %s\n",name);return -1;}

 strcpy(name,argv[1]); /*input file name*/
 strcpy(fname,argv[2]); /*output file name*/
 nsub2=atoi(argv[3]);  /*sqrt subnets per node*/
 nproc2=atoi(argv[4]); /*sqrt nodes*/
 nsub=nsub2*nsub2;    /* number of subnets per node */
 nproc=nproc2*nproc2; /* number of nodes */
 
 fscanf(inputfile,"%d",&torus_flag);
 fscanf(inputfile,"%d",&lrodz);

 cfg=fopen(strcat(name,".cfg"),"r");
 if(cfg==NULL)
    {fprintf(stderr,"Can't open configuration file %s\n",name);return 2;}

  set_flags();
  fscanf(cfg,"%d",&dummy);
  fscanf(cfg,"%d",&dummy);
  fscanf(cfg,"%s",&strdummy);
  fscanf(cfg,"%d",&dummy);
  fscanf(cfg,"%s",&strdummy);
  fscanf(cfg,"%d",&dummy);
  fscanf(cfg,"%d,%d",&dummy,&histoint);
  /* printf("dummy=%d,histoint=%d\n",dummy,histoint); */
 
  fclose(cfg);
 
  flag=(int**)calloc(nsub,sizeof(int*));
  for(i=0;i<nsub;i++)
    flag[i]=(int*)calloc(nproc,sizeof(int));  /* Subnet# X Node# flag array */ 

  flagv=(int**)calloc(nsub,sizeof(int*));     /* Subnet# X Node# flag array */
  for(i=0;i<nsub;i++)
    flagv[i]=(int*)calloc(nproc,sizeof(int));
 
                                       /* table for rodzajow neurons  */
 N=(int**)malloc(lrodz*sizeof(int*));        /* table allocation for size of networks for each class  */
  for(i=0;i<lrodz;i++)
    N[i]=(int*)malloc(3*sizeof(int));
 k=(int**)malloc(lrodz*3*sizeof(int));        /* table allocation for network scaling factor for each class  */
  for(i=0;i<lrodz;i++)
    k[i]=(int*)malloc(3*sizeof(int)); 

 o=(int ***)malloc(lrodz*sizeof(int **));  /* allocating neighborhood tables  */
   for(i=0;i<lrodz;i++)
     o[i]=(int **)malloc(lrodz*sizeof(int *));
    for(i=0;i<lrodz;i++)  
     for(j=0;j<lrodz;j++)
        o[i][j]=(int *)malloc(3*sizeof(int));

 s=(int **)malloc(lrodz*sizeof(int*));           /* synaptic table allocation */
   for(i=0;i<lrodz;i++)
    s[i]=(int*)malloc(lrodz*sizeof(int));

   w=(int **)malloc(lrodz*sizeof(int*));           /* allocation for synaptic weight storage */
   for(i=0;i<lrodz;i++)
    w[i]=(int*)malloc(lrodz*sizeof(int));

 sdw=(int **)malloc(lrodz*sizeof(int*));           /* allocation for synaptic weight dispersion  */
   for(i=0;i<lrodz;i++)
    sdw[i]=(int*)malloc(lrodz*sizeof(int));

 d=(int **)malloc(lrodz*sizeof(int*));           /* allocation for synaptic delay */
   for(i=0;i<lrodz;i++)
    d[i]=(int*)malloc(lrodz*sizeof(int));

 sdd=(int **)malloc(lrodz*sizeof(int*));           /* allocation for synaptic delay dispersion */
   for(i=0;i<lrodz;i++)
    sdd[i]=(int*)malloc(lrodz*sizeof(int));


/* input parameters from files to fill arrays */

strcpy(buff,fname); /* copying output file name to 'buff'*/
nodefile = fopen(strcat(buff,".node"),"wt"); /* file used by jhist program for linking histograms */
fprintf(nodefile,"%d\n",lrodz);

for(i=0;i<lrodz;i++)
  {
   fscanf(inputfile,"%d,",&nr);
   fprintf(nodefile,"%d\n",nr);
     N[i][0]=N[i][1]=nr; /* square mapping */
     N[i][2]=1;
     k[i][0]=k[i][1]=k[i][2]=1;
  }
 fprintf(nodefile,"%d\n",nsub2);
 fprintf(nodefile,"%d\n",nproc2);
 fclose(nodefile);
 fscanf(inputfile,"%d",&nr); /* reading summation check on total number of neurons per subnetwork */

/* setting neighborhood and other input parameters  */

 for(i=0;i<lrodz;i++)
  for(j=0;j<lrodz;j++)
   { 
     fscanf(inputfile,"%d,%d,%d",&o[j][i][0],&o[j][i][1],&o[j][i][2]); /* connection range for given neuron */
     fscanf(inputfile,"%d",&s[j][i]);  /* number of links from given neuron */
     fscanf(inputfile,"%d",&w[j][i]);  /* synaptic weights for given neuron */
     fscanf(inputfile,"%d",&sdw[j][i]);/* weight dispersion */
     fscanf(inputfile,"%d",&d[j][i]);  /* synaptic delay */
     fscanf(inputfile,"%d",&sdd[j][i]);/* delay dispersion */
   }
 
 fclose(inputfile);

/* Setting up noise flags in array */

 iproc=(nproc2/2)*(nproc2+1); /* works for even case */


if(nproc%2)
   isub=0;
 else
   isub=(nsub2/2)*(nsub2+1);


/* Old Central Noise Source Routine */
//if(isub=0) /* central noise source for oddXodd node array including case of one processor */
//   {
//     isub=(nsub2/2)*(nsub2+1);
//     iproc = (nproc-1)/2;
//     flag[isub-nsub2-1][iproc]=flag[isub-nsub2][iproc]=flag[isub-1][iproc]=flag[isub][iproc]=4;  // 4 neurons marked with noise flag
//     printf("4 neurons with NOISE in subnets %d,%d,%d,%d in node %d\n",isub-nsub2-1,isub-nsub2,isub-1,isub,iproc);
//   }
// else  /*sets central four neurons with noise flags in evenXeven case*/
//   {
//   flag[nsub-1][iproc-nproc2-1]=flag[nsub-nsub2][iproc-nproc2]=flag[nsub2-1][iproc-1]=flag[0][iproc]=4; 
//    printf("4 neurons with NOISE in subnet/node %d/%d,%d/%d,%d/%d,%d/%d\n",
//   	   nsub-1,iproc-nproc2-1,nsub-nsub2,iproc-nproc2,nsub2-1,iproc-1,0,iproc);
   //for(i=0;i<nsub;i++)
   //  for(j=0;j<nproc;j++)
   //    flag[i][j]=4;
//   }


/* Setting up vout flags in array for use with SVE program */
   if(nsub2%2)
     vsub=(nsub2*(nsub2-1))/2+(nsub2-1)/2; /* for odd, provides index of middle subnetwork */
   else
   vsub=nsub2*nsub2/2-nsub2/2;
 
   flagv[vsub][0]=flagv[vsub][nproc2-1]=flagv[vsub][nproc2*(nproc2-1)]=flagv[vsub][nproc-1]=nr; /* nr is the total number of neurons/subnet */  

/*

OFFSET CALCULATIONS TO FILL TABLES

*/

 NN=(int*)malloc(lrodz*sizeof(long int));   /* table with numbers of neurons within each class */
 for(i=0;i<lrodz;i++)
  NN[i]=(N[i][0]?N[i][0]:1)*(N[i][1]?N[i][1]:1)*(N[i][2]?N[i][2]:1);
 
 NOF=(int*)malloc((lrodz+1)*sizeof(int));  /* table with offsets */
  for(i=1,NOF[0]=0;i<=lrodz;i++)
    NOF[i]=NOF[i-1]+NN[i-1];

 NB=(int*)calloc(lrodz,sizeof(int));     /* base no of neurons in subnets */
  if (NB==NULL){fprintf(stderr,"Insufficient memory for  NB\n");exit(1);}

 

/* maximal neighborhood calculation */

 maxot=0;
 if(torus_flag!=2)
 for(i=0;i<lrodz;i++)
  for(j=0;j<lrodz;j++)
   {
     tmpot=1;
     for(l=0;l<3;l++)
     tmpot*=2*o[i][j][l]+1;
      if(tmpot > maxot)
       maxot=tmpot;
   }
  else
   for(i=0;i<lrodz;i++)
    {
    if(maxot<NN[i])
       maxot=NN[i];
    } 
 tab=(int*)malloc(maxot*sizeof(int));

 Nneu=0;
 for(i=0;i<lrodz;i++)
     Nneu+=NN[i];   /* Calculate the total number of neurons in the subnetwork */

 if(Nneu!=nr) {fprintf(stderr,"\nerror: check lines 2 & 3 in input file\n");
   return 2;
   }else
     Nneu*=nsub; /* Total number of neurons of all classes per node */

 if(Nneu>MAX_NEURON)
   {fprintf(stderr,"No of Neurons per node : %d > %d\n",Nneu,MAX_NEURON);
   return 2;
   }
 //printf("Nneu=%i\n",Nneu);

 neu=(struct NEURON *)calloc(Nneu,sizeof(struct NEURON));
   if (neu==NULL){fprintf(stderr,"Insufficient memory for neu");exit(1);}

   for(iproc=0;iproc<nproc;iproc++) /* nproc is number of nodes, big loop over nodes */
   {
     ulong nsum=0;

     sprintf(name,"%s.syncomp.%d",fname,iproc);
     if(!(synfile=fopen(name,"w")))
	  {fprintf(stderr,"Can't open %s\n",fname);return 2;}
     sprintf(name,"%s.neutmp.%d",fname,iproc);
     if(!(neufile=fopen(name,"w")))
          {fprintf(stderr,"Can't open %s\n",fname);return 2;}
     for(i=0;i<lrodz;i++)
	 NB[i]=NOF[i]*(nsub-1);
     ns=0;                           /* number of synapses, talleyed below when writing synapse file  */

     sprintf(name,"tmpsyn.%d",iproc);
     if(!(outputfile=fopen(name,"wb")))
          {fprintf(stderr,"Can't open %s\n",name);return 2;}
     //outputfile=fopen("syn.tmp","wb");
     //if(outputfile==NULL){fprintf(stderr,"\nError in opening output file");exit(1);}

     for(isub=0;isub<nsub;isub++)
        {
	 
          //for(i=0;i<lrodz;i++)
          //printf("\nNN[%d]=%d, NOF[%d]=%d",i,NN[i],i,NOF[i]); /* checks */
          //printf("\nmaxot=%d,Nneu=%d\n",maxot,Nneu);
     
          /*

          3.Neighborhood calculations

          */
 

          /*randomize();*/
          nr=0;                           /* neighborhood counting */

          for(ir=0;ir<lrodz;ir++)  /* loop over neuron class */
             for (ix=0,xyz[0]=0;ix<N[ir][0];ix++,xyz[0]+=k[ir][0])     /* x-loop */
                for(iy=0,xyz[1]=0;iy<N[ir][1];iy++,xyz[1]+=k[ir][1])    /* y-loop */
                   for(iz=0,xyz[2]=0;iz<N[ir][2];iz++,xyz[2]+=k[ir][2],nr++)  /* z-loop  */
                      for(jp=0;jp<lrodz;jp++)     /* loop on classes of neurons */
                      {
                              /*     wygeneruj tablice z numerami neuronow rodzaju jp 
                                     w otoczeniu neuronu nr (czyli neuronu rodzaju ir,
                                     w polozeniu ix,iy,iz) np.*/
			  switch(torus_flag) /* tabulating number of neurons available for connection */
                              {
                                  case 0:
				     no=otoczenie(nr,ir,xyz,jp);  /* no torus boundary condition */
                                     break;
                                  case 1:
                                     no=otoczenie2(nr,ir,xyz,jp); /* torus boundary condition */
                                     break;
                                  case 2:                         /* case 2 choosing from whole array  */
                                     for(i=NOF[jp],no=0;i<NN[jp]+NOF[jp];i++)
                                     if(i!=nr)
                                     tab[no++]=i;
                                     break;
                                  default:
                                     fprintf(stderr,"wrong option for torus flag %d\n",torus_flag);
                                     exit(0);
                              }


			  if(!no)continue;
			  l=(no<s[ir][jp])?no:s[ir][jp];
			  /* random connection formation */
			  for(i=0;i<l;i++) 
			  {
			      /* random connection formation */
			      nl=losuj(tab,no); /*printf("%d,%d,%d,%d\n",ir,jp, nr,nl);*/
			      /* filling of tables for delay/synaptic weights */
			      ns+=zapisz(nl,nr,ir,jp,outputfile); /*writing output to file 0.neu, 0.syn  */    

			      /* printf("\npostindex=%d, postclass=%d, preindex=%d, preclass=%d",nr,ir+1,nl,jp+1); */ /* print check variables */
                 
			  } /* end of drawing */
                      } /* jp */
	              for(i=0;i<lrodz;i++)
	              NB[i]+=NN[i];

        } // end of isub

  fclose(outputfile);
 

  sprintf(name,"tmpsyn.%d",iproc);
  if(!(outputfile=fopen(name,"r")))
       {fprintf(stderr,"Can't open %s\n",name);return 2;}
  //outputfile=fopen("syn.tmp","r");

  syn_tmp=(struct SYN_TMP *)calloc(ns,sizeof(struct SYN_TMP)); /* tablica wszystkch syn_tmp */
   if (syn_tmp==NULL){fprintf(stderr,"Can't allocate temporary synapse memory space\n");exit(1);}
  fread(syn_tmp,sizeof(struct SYN_TMP),(size_t)ns,outputfile); /* czytaj syn_tmp z syn.tmp */   
  fclose(outputfile);
  //printf("\nns=%d iproc=%d\n",ns,iproc);/* print check variables */
  //unlink(name);
  //unlink("syn.tmp");  /* syn.tmp nie jest dalej potrzebny */

  qsort( (void *)syn_tmp,(size_t)ns,sizeof(struct SYN_TMP),sort_neu_num ); /* sorting synapses by presynaptic neuron index for simulator use
									    sort_neu_num defined below*/ 
  
  neu_pre=calloc(Nneu,sizeof(int));                 /* table of neurons with number of synapses (presynaptic neurons)*/
  is=0;syn_n=0;
  for(i=0;i<Nneu;i++)
   {
    neu_pre[i]=0;
    while(i==(syn_tmp+syn_n)->neuron_pre){
    syn_n++;
    neu_pre[i]++; }
    qsort((void *)(syn_tmp+is),(size_t)neu_pre[i],sizeof(struct SYN_TMP), sort_delay );   /* sort by delay */
    is=syn_n;
   } 



  for(is=0;is<ns;is++)
   {
     //     int kind=get_kind((syn_tmp+is)->neuron_post);
    syn.neuron=(syn_tmp+is)->neuron_post;
    syn.weight=(syn_tmp+is)->weight;
    syn.delay=(syn_tmp+is)->delay;
    //printf("=%d,iproc=%d,i=%d,j=%d\n",isub,iproc,i,j);
    syn.preclass_num=(syn_tmp+is)->preclass_num;
    fwrite(&syn,sizeof(struct SYNAPSE),1,synfile);
   }

  nsum+=ns;
  free(syn_tmp);
 
  /* teraz neurony */
 
  for(i=0;i<Nneu;i++){
     (neu+i)->syn_num=neu_pre[i];   /* preparing to write neuron files, structure NEURONS */
     set_neu(neu+i);}        /* setting initial voltage and other values */

  

  for(i=0;i<lrodz;i++)
     for(j=0;j<NN[i]*nsub;j++)
        (neu+NOF[i]*nsub+j)->param=i;      /* assigning type of neuron */

  /* noise flags */
  //for(isub=0;isub<nsub;isub++)
  //   { l=isub*NN[0];
  //      if(flag[isub][iproc]){      
  //         for(i=0;i<flag[isub][iproc];i++)
  //         {
 	      //j=randomm(NN[0])+l;
	      //j=l+N[0][0]/2*(N[0][0]+1)+i;
  //	      j=l+i; 
              //j=isub*N[0][0]*N[0][0]+i;
  //            printf("Noise flags ** isub=%d,iproc=%d,i=%d,j=%d\n",isub,iproc,i,j);
  //	      (neu+j)->flags|=NOISE;  /* only 1st kind with NOISE */
  //         } 
  //      }    /* end if */
  //   }


  /* noise out flags */
  gumba=0;
  for(isub=0;isub<nsub;isub++)
  { 

           l=0;
           for(n=0;n<lrodz;n++)
           {
              for(i=0;i<NN[n];i++)  
              {
                 j=l+isub*NN[n]+i;
                 //printf("Noise flags ** isub=%d,iproc=%d,i=%d,j=%d\n",isub,iproc,i,j);
		 testrandnoise=(double)random()/((double)RAND_MAX+1);
		 // printf("testrandnoise=%f\n",testrandnoise);
		 if( (testrandnoise<0.01) && ( (n==0) || (n==2) ) ){
                      (neu+j)->flags|=NOISE;
		 }

              } 
              l+=nsub*NN[n];
            }
         
  }

  /* vout flags */
  for(isub=0;isub<nsub;isub++)
     { 

        //l=isub*NN[0];
        //m=nsub*NN[0]+isub*NN[1];  

        if(flagv[isub][iproc])
        {
           l=0;
           for(n=0;n<lrodz;n++)
           {
              for(i=0;i<NN[n];i++)  
              {
                 j=l+isub*NN[n]+i;
                 printf("Voltage out flags ** vsub=%d,iproc=%d,i=%d,j=%d,class=%d\n",isub,iproc,i,j,n);
                 (neu+j)->flags|=VOUT;  /* only 1st kind VOUT */
                 //(neu+j)->flags|=IEXT; 
              } 
              l+=nsub*NN[n];
            }
         }

         /*
         for(i=0;i<flagv[isub][iproc];i++)
         {
            j=l+i;
            printf("vsub=%d,iproc=%d,i=%d,j=%d\n",isub,iproc,i,j);
            (neu+j)->flags|=VOUT;  // only 1st kind VOUT 
         }

         for(i=0;i<sqrt(flagv[isub][iproc]);i++) // inhibitorty 
         {
            j=m+i;
            printf("vsub=%d,iproc=%d,i=%d,j=%d\n",isub,iproc,i,j);
            (neu+j)->flags|=VOUT;  // 2nd kind VOUT 
         }
         */
     }
 
  fwrite(neu,sizeof(struct NEURON),Nneu,neufile);

 
  fclose(synfile);
  fclose(neufile);
  
  sprintf(name,"%s.cfg.%d",fname,iproc);
  cfg=fopen(name,"w");   /* zbior *.cfg */
  //fprintf(outputfile,"%s\n","ONSET");         /* creating .cfg file */
  fprintf(cfg,"%s\n",flags);                  /* */
  fprintf(cfg,"%s\n",fname);                  /* */
  fprintf(cfg,"%d\n",lrodz);
  fprintf(cfg,"%s\n",fname);                  /* */
  fprintf(cfg,"%lu\n",nsum);             /* */
  fprintf(cfg,"%s\n",fname);                  /* */   
  fprintf(cfg,"%d\n",Nneu);            /* */
  fprintf(cfg,"%d,%d\n",Nneu,histoint);      /* histogram parameters-length of bin  */ 
  fprintf(cfg,"0\n");                /* out_no*/
  fprintf(cfg,"0.00003068");            /* the lambda actually used by simulator, was 3.068X10^-4 */
  fclose(cfg); 
 } // end iproc


    /* copy par file */

  strcpy(name,argv[1]);
  par=fopen(strcat(name,".par"),"r");
  if(par==NULL)
    {fprintf(stderr,"Can't open %s\n",name);return 2;}
 
  for(i=0;i<nproc;i++)
    {
      int ni=0;
      sprintf(name,"%s.par.%d",argv[2],i);
      opar=fopen(name,"w");
      if(opar==NULL)
	{fprintf(stderr,"Can't open %s\n",name);return 2;}
      rewind(par);
      while(ni<lrodz && fread(&p,sizeof(struct PARAMETERS),1,par)==1 && 
	    fread(&delays[ni][0],sizeof( SPIKE_INTERVAL ),lrodz,par)==lrodz)
	{
	     
	  fwrite(&p,sizeof(struct PARAMETERS),1,opar);
	  fwrite(&delays[ni][0],sizeof( SPIKE_INTERVAL ),lrodz,opar);
	  ni++;
	}
      fclose(opar);
    }
  fclose(par);

  /* copy stim file */
  if(stim_f)
     {
        strcpy(name,argv[1]);
        par=fopen(strcat(name,".stim"),"r");
        if(par==NULL)
	   {fprintf(stderr,"Can't open %s\n",name);return 2;}
    
        for(i=0;i<nproc;i++)
	{
	   rewind(par);
	   sprintf(name,"%s.stim.%d",argv[2],i);
	   opar=fopen(name,"w");
	   if(opar==NULL)
	      {fprintf(stderr,"Can't open %s\n",name);return 2;}
	   while(!feof(opar))
	      {  char buf[100];
	         if(fscanf(par,"%[^\n]\n",buf)!=1)break;
	         fprintf(opar,"%s\n",buf);
	      }
	    fclose(opar);
	      
	 }
         fclose(par);
      }
  

printf("\n");
return 0;
} /*koniec main()*/


/*
sort_delay  - needed in qsort
*/

int sort_delay(const void * a, const void * b )
{
char delay1[1],delay2[1];
struct SYN_TMP *syn1=NULL, *syn2=NULL;
syn1=(struct SYN_TMP *)a;
syn2=(struct SYN_TMP *)b;
*delay1=syn1->delay;
*delay2=syn2->delay;

return( strcmp(delay1,delay2) );
}


/*
sort_neu_num - needed in qsort() 
*/

int sort_neu_num(const void * a, const void * b)
{
NEURON_INDEX neu1,neu2;
struct SYN_TMP *syn_tmp1=NULL,*syn_tmp2=NULL;
syn_tmp1=(struct SYN_TMP *)a;
syn_tmp2=(struct SYN_TMP *)b;
neu1=syn_tmp1->neuron_pre;
neu2=syn_tmp2->neuron_pre;
return(neu1-neu2);
}

/*
  
4. Generation of neighborhood table:

*/

int rozrzut(int sd)
{
 int los;
 los=randomm(100);
  if(los < 50 )
   return( -1*randomm(sd) );
  else
   return(    randomm(sd) );
}


int losuj(int *tablica,int index) /* random connections within identified regions of allowed connections */
{
int i,j;
do
{
i=randomm(index);
}while(tablica[i]<0);
j=tablica[i];
tablica[i]=-1;
return j;
}

int zapisz(int inr, int inl, int i, int j, FILE * outputfile)
{
struct SYN_TMP syn_tmp; 
extern  int rozrzut(int);

syn_tmp.neuron_post=inl+NB[i];
syn_tmp.neuron_pre=inr+NB[j]; 
syn_tmp.preclass_num=j;      
syn_tmp.weight=w[i][j]+rozrzut(sdw[i][j]);
   /* obliczamy najpierw wspolrzedne i[3] neuronu na siatce jp, ktory jest najblizej x,y,z */
syn_tmp.delay =d[i][j]+rozrzut(sdd[i][j]);
//printf("pre= %d,post= %d,preclass= %d, weight= %d,delay= %d\n",syn_tmp.neuron_pre,syn_tmp.neuron_post,syn_tmp.preclass_num,syn_tmp.weight,syn_tmp.delay);
return(fwrite(&syn_tmp,sizeof(struct SYN_TMP),1,outputfile));
}


int otoczenie(int nr ,int ir,int xyz[3],int jp)
{
   int nx,ny,no,i[3],ix,iy,iz,mini[3],maxi[3], *ot;
   int j;

   /* tab must be declared earlier with sufficient size */

   ot = o[ir][jp];
   
   /*   The next line calculates the coordinates of neurons from one class to the others */  

  no=0;
  for(j=0;j<3;j++)
         i[j]=xyz[j]/k[jp][j];  

  /* No torus boundary condition */

  for(j=0;j<3;j++)
  {
     mini[j]=i[j]-ot[j];
     maxi[j]=i[j]+ot[j]+1;
     if(mini[j]<0)mini[j]=0;
     if(maxi[j]>N[jp][j])maxi[j]=N[jp][j];
  }
   
   for(ix=mini[0];ix<maxi[0];ix++)
      { nx=ix*N[jp][1]*N[jp][2]+NOF[jp];
      for(iy=mini[1];iy<maxi[1];iy++)
         { ny=nx+iy*N[jp][2];
         for(iz=mini[2];iz<maxi[2];iz++)
         {
            if((j=ny+iz)!=nr)
            tab[no++]=j;   
         }
      }
   }

return no;
}


int otoczenie2(int nr, int ir,int xyz[3],int jp)
{
   int nx,ny,no,i[3],ix,iy,iz,mini[3],maxi[3], *ot;
   int j,ii,NX,NY,NZ;
   /* Torus boundary condition */

   ot = o[ir][jp];
   /* calculating location of neurons used to create connections  x,y,z */
  
  
  no=0;
  for(j=0;j<3;j++)
         i[j]=xyz[j]/k[jp][j];  

  /* teraz otoczenie
     wersja z zawijaniem na torusie */

      
 for(j=0;j<3;j++)
 {
   mini[j]=i[j]-ot[j];
   maxi[j]=i[j]+ot[j]+1;
 }
   NX=N[jp][0];
   NY=N[jp][1];
   NZ=N[jp][2];
 
   for(ix=mini[0];ix<maxi[0];ix++)
   {  ii=(ix+NX)%NX;
      nx=ii*NZ*NY+NOF[jp];
      for(iy=mini[1];iy<maxi[1];iy++)
      {  ii=(iy+NY)%NY;
         ny=nx+ii*NZ;
         for(iz=mini[2];iz<maxi[2];iz++)
         {
             if((j=ny+(iz+NZ)%NZ)!=nr)
             tab[no++]=j;
         }
      }
    }
return no;
} 


void set_neu(struct NEURON * neu)
{
   int j,kk;
   switch(neu->param)
   {
      default:
      neu->flags = HISTO;
      neu->polarity = 0x0;

      for (kk=0;kk<NEXCCLASS;kk++){	

           neu->Ve_o[kk]=0;
           neu->Ve_d[kk]=0;

      }


      for (kk=0;kk<NINHCLASS;kk++){
	
           neu->Vi_o[kk]=0;
           neu->Vi_d[kk]=0;

      }


      neu->interval=0;
      neu->sum_interval=0;
      for(j=0;j<SPIKE_LIST_LENGTH;j++)
         neu->spikes[j]=0;
      neu->first=EMPTY_LIST;
      neu->last=0x1;
      neu->V=V_0;
      neu->W=W_0;
#ifdef CALCIUM
      for(j=1;j<7;j++)
      {     
         neu->C[j]=0.;
         neu->U[j]=0.;
      } 
//      neu->C=C_0;
      neu->X=X_0;
      neu->B=B_0;
#endif
   }
 
}

unsigned seed_save()
{
 unsigned seed;
 seed = (unsigned)time(NULL);
 srandom(seed);
return seed;
}