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]
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 V1 L6 pyramidal corticothalamic cell; Neocortex V1 L2/6 pyramidal intratelencephalic cell; Neocortex V1 interneuron basket PV 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 V1 L6 pyramidal corticothalamic cell; Neocortex V1 L2/6 pyramidal intratelencephalic cell; Neocortex V1 interneuron basket PV cell; AMPA; Gaba; I A; I K; I K,leak; I K,Ca; I Sodium; I Calcium;
/* Program do generacji polaczen synaptycznych miedzy subnets
   Na wyjsciu ma produkowac plik z polaczeniami z zapisana para neuronow waga i
   delay. Nastepny program (lub druga czesc tego samego) moze wykorzystac to
   do generacji struktury neuronow i synaps we wlasciwym formacie.

   01/29/03 P. Kudela, losuj2 prawdopodobienstwo polaczenia zanika exp.  
                         a w zapisz2 delay wzrasta liniowo  
   05/29/03 P. Kudela, now multi-class supported 
*/

#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;               /* number of neuron classes */
int *NN;                 /* NN[lrodz] total number of neurons within each class */
int *NOF;                /* NOF[lrodz] table of index offsets */ 

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 */
double **Ap,**lambda;/* Ap[lrodz][lrodz] lambda[][] amplitude and spatial decay rate of connection probability  */ 

int *tab;         /* temporary neighborhood table tab[maxot]*/ 
double *dist, *costheta, *sintheta;      /* table of distances between candidate neurons and costheta, sintheta, updated 08/23/05 */
int *postx, *posty; /* Array of postsynaptic neuron coordinates for accepted target */
double *zclass;  /* Array of class z-values within cortical model */
int *horaxonflagpos,*veraxonflagpos, *axonflag; /* Array of flags for stimulated axon portions */
int *horaxonflagneg,*veraxonflagneg;

/* Neuron Class            Type               Firing Pattern
        1           Layer II/III Pyramid             RS
        2           Layer IV Spiny                   RS
        3           Layer V Pyramid                  IB
        4           Layer VI Pyramid                 IB
        5           Basket Layers II-VI              FS
        6           Double Bouquet II-VI             FS
        7           Chandelier                       FS     */

#define P23 1
#define P4  2
#define P5  3
#define P6  4
#define BASKET 5
#define BOUQUET 6
#define CHANDELIER 7
#define BASKETANG 0.0 //1.5708  /* for offset costheta distr for basket cell axons */
#define PIIIANG 0.0   /* for offset costheta distr for PIII cell axons */
#define PVANG 0.0 /*for offset costheta distr for PV cell axons */
#define PVIANG 0.0    /* for offset costheta distr for PVI cell axons */

#define SPATIALSCALE 25 /* center to center column spacing in microns */

#define BASKETLAMBDA 50 /* 1/e decay for basket cell axon connections */
#define BASKETNEAR 50 /* distance in microns within which basket conn'ions are isotropic */

#define PIIINEARRAD 300   /* Near connection radius for PIII neurons */
#define STRIPWIDTH  500   /* 500 micron strip width for domain simulations */
#define PIIIFARAVE 2000   /* Ave dist to long range connection spot */
#define PIIISPOTRAD 200   /* Long range spot radius */

#define BOUNDARYLAYER 300 /* Boundary layer distance, to reduce connections */

#define ISLENGTH    30    /* One half initial segment length (microns) */

/* Definitions for E-field effects */

#define SEG 500  /* spacing btwn nodes of Ranvier for 5.7 micron myelinated axon */
//#define VOLT .2 /* stimulation voltage applied to metal electrode in voltage controlled mode */
#define I0 .010 /* Stimulation current in current controlled mode in Amps */
#define SIGMA 0.3 /* Tissue conductivity in A/Vm */
#define THRESHD2V 2 /* Threshold for activation fct (mV), second difference of voltage along fiber element */ 
#define A 1000 /* Electrode disk radius (microns) */
#define STIMPOS 2340 /* Electrode Z-value (microns), taken with z=0 at the top of the white matter */
#define DISKX 32 /* Center of electrode disk x-value (lattice constants) */
#define DISKY 32 /* Center of electrode disk y-value (lattice constants) */ 
#define PI 3.14159

/* Cortical layer z-values */
#define ZVIB 162 /* VIbeta */
#define ZVIA 488 /* VIalpha */
#define ZV   1130 /* V */
#define ZIVCBETA 1230 /* IVCbeta */
#define ZIVCALPHA 1230 /* IVCalpha */
#define ZIVB 1230 /* IVB */
#define ZIVA 1230 /* IVA */
#define ZIII 1540 /* III */
#define ZII  1940 /* II */
#define ZI   2185 /* I */

#define Gi   7.6e-07 /* Intracompartment conductivity in inverse Ohms */ 

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

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

int main(int argc, char  *argv[])
{
 int i,j,l,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[10];
 FILE *inputfile,*outputfile,*graphfile,*axflagfile;
 struct SYN_TMP *syn_tmp;
 struct SYNAPSE syn;
 struct NEURON *neu;
 int otoczenie(int,int,int*,int);
 int otoczenie2(int,int,int*,int);
 int losuj(int*,int);
 int losuj2(int*,double*,double*,double*,int,double,double,int,int,int*);
 extern int zapisz(int,int,int,int,FILE *);
 extern int zapisz2(int,int,int,int,double,FILE *,FILE *,FILE *,int);
 extern int sort_neu_num(const void *, const void *);
 int nsub;
 int pre[3];
 double horzmid,horxmid,horymid,horrhomid,horVmid,horrhopre,horVpre,horrhopost,horVpost;
 double horseglength,hordelsquaredV;
 double verzmid,verxmid,verymid,verrhomid,verVmid,verrhopre,verVpre,verrhopost,verVpost;
 double verseglength,verdelsquaredV,aradius,coeff;
 int synapsesum;
 double currentinjsum, avecurrentinj;
 int nps,s_j,s_i,critdex;

/*

1.Reading and checking input file information

*/

/* printf("Made it to file checking\n"); */

 if(argc<3)
   {
     fprintf(stderr,"USAGE:\n\tmk_link <file_name> <sqrt(no of subnets in all)> [<seed>]\n");
     return 1;
   }
 if(argc>3)
   srandom((unsigned)atoi(argv[3]));
  else
   printf("seed %u\n",seed_save()); /* seed for srand */
 srand48(127L);
 //inputfile = fopen(strcpy(name,argv[1]),"r");
 inputfile = fopen(strcat(strcpy(name,argv[1]),".net"),"r");
 if(inputfile==NULL){fprintf(stderr, "Network definition file not present %s\n",name);return -1;}

 fscanf(inputfile,"%d",&torus_flag);
 fscanf(inputfile,"%d",&lrodz);

 nsub=atoi(argv[2]);
                                       /* Memory allocations  */
 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 weights */
   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));

  Ap=(double **)malloc(lrodz*sizeof(double));      /* used to compute probability of connection */
   for(i=0;i<lrodz;i++)
    Ap[i]=(double*)malloc(lrodz*sizeof(double));

  lambda=(double **)malloc(lrodz*sizeof(double));  /* used to commpute probability of connection */
   for(i=0;i<lrodz;i++)
    lambda[i]=(double*)malloc(lrodz*sizeof(double));

  zclass=(double*)malloc(lrodz*sizeof(double));

/* Load cortical layer z-information for E-field calc, array index by class number */
   zclass[0]= ZII;
   zclass[1]= ZIVB;   
   zclass[2]= ZV;  
   zclass[3]= ZVIA;
   zclass[4]= ZII;  
   zclass[5]= ZII;  
   zclass[6]= ZIVB;
 
/* initialize mapping and scaling arrays */

for(i=0;i<lrodz;i++)
{
   N[i][0]=N[i][1]=nsub; //square mapping
   N[i][2]=1;
   k[i][0]=k[i][1]=k[i][2]=1;
}


/* input parameters from file to fill arrays */


 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]);
     fscanf(inputfile,"%d",&s[j][i]);
     fscanf(inputfile,"%d",&w[j][i]);
     fscanf(inputfile,"%d",&sdw[j][i]);
     fscanf(inputfile,"%d",&d[j][i]);
     fscanf(inputfile,"%d",&sdd[j][i]);
     fscanf(inputfile,"%lf",&Ap[j][i]);    
     fscanf(inputfile,"%lf",&lambda[j][i]); 
   }
 fclose(inputfile);
 

/*

2. Offset calculations and table filling

*/

 /* printf("Made it to offset calculations\n"); */

 NN=(int*)malloc(lrodz*sizeof(long int));   /* table with number of subnetworks containing 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*sizeof(int)); /* table with offsets */
  for(i=1,NOF[0]=0;i<lrodz;i++)
    NOF[i]=NOF[i-1]+NN[i-1];


/* maximal neigborhood calculation */

 maxot=0;
 if(torus_flag!=2) /* i.e. not choosing from whole array */
 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));
 dist=(double*)malloc(maxot*sizeof(double));
 costheta=(double*)malloc(maxot*sizeof(double));
 sintheta=(double*)malloc(maxot*sizeof(double));
 postx=(int*)malloc(maxot*sizeof(int));
 posty=(int*)malloc(maxot*sizeof(int));
 horaxonflagpos=(int*)malloc(maxot*sizeof(int));
 veraxonflagpos=(int*)malloc(maxot*sizeof(int));
 horaxonflagneg=(int*)malloc(maxot*sizeof(int));
 veraxonflagneg=(int*)malloc(maxot*sizeof(int));
 axonflag=(int*)malloc(maxot*sizeof(int));

 Nneu=0;
 for(i=0;i<lrodz;i++)
  Nneu+=NN[i];   /* calculate subnetworks X # of classes */

 for(i=0;i<lrodz;i++)
  printf("\nNN[%d]=%d, NOF[%d]=%d",i,NN[i],i,NOF[i]); /* co by sprawdzic */
  printf("\nmaxot=%d,Nneu=%d\n",maxot,Nneu);
 
/*

3.Setting up for output of links

*/
 
 outputfile=fopen("syn.tmp","wb");
 graphfile=fopen("syn.graph","wb");
 axflagfile=fopen("syn.axflag","wb");

 if(outputfile==NULL){fprintf(stderr,"\nError in opening synapse output file");exit(1);}
 if(graphfile==NULL){fprintf(stderr,"\nError in opening synapse output file");exit(1);}
 if(axflagfile==NULL){fprintf(stderr,"\nError in opening axflagfile file");exit(1);}


    /*randomize();*/
     nr=0;                           /* subnetwork "neuron" number */
     ns=0;                           /* synapse index */
     currentinjsum=0;                /* in microamps */
     synapsesum=0;
     for(ir=0;ir<lrodz;ir++)         /* loop on classes of pre-neurons */
     {
	 //printf("Working on Preclass %d\n",ir+1);
        for (ix=0,xyz[0]=0;ix<N[ir][0];ix++,xyz[0]+=k[ir][0])     /* x-loop */
	{
	    //printf("ix= %d\n",ix+1);
           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 post-neurons */
                  {
                     /*     generate tables with numbers of neurons of each class jp,in the neighborhood of neuron nr (neuron from class ir
			    connects to jp (used to be the other way around) */
                     switch(torus_flag)
                     {
                        case 0:                           /* no torus boundary condition */
                           no=otoczenie(nr,ir,xyz,jp);
                           break;
                        case 1:                           /* torus boundary */
                           no=otoczenie2(nr,ir,xyz,jp);  
                           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[jp][ir])?no:s[jp][ir];
                     /* if no is greater than the number of allowed synapses, then set l=# of allowed synapses  */
		     //printf("Allowed/required conn comp no=%d,preclass=%d,postclass=%d,s[jp][ir]=%d\n",no,ir+1,jp+1,s[jp][ir]);/* network linkage checks */


		     /* Boundary effect correction, reducing allowed connections for neurons in specific column number ranges, and row number ranges */

		     nps=nsub; //sqrt(subnets/node)*sqrt(nodes)
                     s_j=(nr)%nps; /* column index of presynaptic subnetwork */
                     s_i=(nr)/nps-ir*nps; /* row index of presynaptic subnetwork */
		     critdex=BOUNDARYLAYER/SPATIALSCALE;

		     //printf("Subnet indices: column=%d row=%d nps=%d critdex=%d nr=%d\n",s_j,s_i,nps,critdex,nr);

		     //printf("Pre-l: l=%d\n",l);

		     if ( (s_j<critdex) && (s_i>=critdex) && (s_i<=((nps-1)-critdex)) ){
		          l=(l*5)/8;
		     }	 

		     if ( (s_j>((nps-1)-critdex)) && (s_i>=critdex) && (s_i<=((nps-1)-critdex)) ){
		          l=(l*5)/8;
		     }	 

		     if ( (s_i<critdex) && (s_j>=critdex) && (s_j<=((nps-1)-critdex)) ){
		          l=(l*5)/8;
		     }	 

		     if ( (s_i>((nps-1)-critdex)) && (s_j>=critdex) && (s_j<=((nps-1)-critdex)) ){
		          l=(l*5)/8;
		     }	 

		     if ( (s_j<critdex) && (s_i<critdex) ){
		          l=(l*3)/8;
		     }	 

		     if ( (s_j<critdex) && (s_i>((nps-1)-critdex)) ){
		          l=(l*3)/8;
		     }	 

		     if ( (s_j>((nps-1)-critdex)) && (s_i<critdex) ){
		          l=(l*3)/8;
		     }	 

		     if ( (s_j>((nps-1)-critdex)) && (s_i>((nps-1)-critdex)) ){
		          l=(l*3)/8;
		     } 

		     //printf("Post-l: l=%d\n",l);

		     /* Form connections and determine stimulation locations */

                     for(i=0;i<l;i++)    /* losuje z calego otoczenia - nowa wersja */
                     {                    /* z losuj2  01/29/03 */ 
                        
                        for(j=0;j<3;j++)
			     pre[j]=xyz[j]/k[jp][j];  

			// printf("no=%d,preclass=%d,postclass=%d,l=%d,i=%d\n",no,ir+1,jp+1,l,i);  /* network linkage checks */
			nl=-1;
			while (nl==-1){
                             nl=losuj2(tab,dist,costheta,sintheta,no,Ap[jp][ir],lambda[jp][ir],jp,ir,pre);
			}
			//printf("i= %d,l= %d,nl= %d\n",i,l,nl);
			if(nl!=-1)
                        {

	                /* Set synapse flag if gradient requirement met */
                        /* Check E field gradient requirement */
			/* Charged balance, bipolar pulse */
   
                        /* Horizontal axon segment */
			aradius=A;
			aradius=aradius/(1000*1000);
			coeff=I0/(4*PI*SIGMA*aradius);
			//printf("\n aradius=%f, coeff=%f",aradius,coeff);


			axonflag[nl]=0;
		
                        horaxonflagpos[nl]=0;
			horaxonflagneg[nl]=0;
			horzmid=zclass[jp];
                        horxmid=(pre[0]+postx[nl])/2;
			horymid=(pre[1]+posty[nl])/2;
			horrhomid=sqrt(pow((DISKX-horxmid),2)+pow((DISKY-horymid),2))*25;
			horVmid=(coeff)*asin(2*A/(sqrt(pow((STIMPOS-horzmid),2)+pow((horrhomid-A),2))+sqrt(pow((STIMPOS-horzmid),2)+pow((horrhomid+A),2))));
			horrhopre=sqrt(pow((DISKX-pre[0]),2)+pow((DISKY-pre[1]),2))*25;
			horVpre=(coeff)*asin(2*A/(sqrt(pow((STIMPOS-zclass[jp]),2)+pow((horrhopre-A),2))+sqrt(pow((STIMPOS-zclass[jp]),2)+pow((horrhopre+A),2))));
			horrhopost=sqrt(pow((DISKX-postx[nl]),2)+pow((DISKY-posty[nl]),2))*25;
			horVpost=(coeff)*asin(2*A/(sqrt(pow((STIMPOS-zclass[jp]),2)+pow((horrhopost-A),2))+sqrt(pow((STIMPOS-zclass[jp]),2)+pow((horrhopost+A),2))));
			horseglength=sqrt(pow((postx[nl]-pre[0]),2)+pow((posty[nl]-pre[1]),2))*25;
			/* hordelsquaredV in mV */ 
			hordelsquaredV = (horVpre+horVpost-2*horVmid)*1000;
                        if (hordelsquaredV > THRESHD2V)
			    horaxonflagpos[nl]=1;
			if (-hordelsquaredV > THRESHD2V)
			    horaxonflagneg[nl]=1;

                        // Vertical axon segment  (Old Version)
                        //veraxonflagpos[nl]=0;
			//veraxonflagneg[nl]=0;
			//verzmid=(zclass[ir]+zclass[jp])/2;
                        //verxmid=pre[0];
			//verymid=pre[1];
			//verrhomid=sqrt(pow((DISKX-verxmid),2)+pow((DISKY-verymid),2))*25;
			//verVmid=(coeff)*asin(2*A/(sqrt(pow((STIMPOS-verzmid),2)+pow((verrhomid-A),2))+sqrt(pow((STIMPOS-verzmid),2)+pow((verrhomid+A),2))));
			//verrhopre=sqrt(pow((DISKX-pre[0]),2)+pow((DISKY-pre[1]),2))*25;
			//verVpre=(coeff)*asin(2*A/(sqrt(pow((STIMPOS-zclass[ir]),2)+pow((verrhopre-A),2))+sqrt(pow((STIMPOS-zclass[ir]),2)+pow((verrhopre+A),2))));
			//verrhopost=sqrt(pow((DISKX-pre[0]),2)+pow((DISKY-pre[1]),2))*25;
			//verVpost=(coeff)*asin(2*A/(sqrt(pow((STIMPOS-zclass[jp]),2)+pow((verrhopost-A),2))+sqrt(pow((STIMPOS-zclass[jp]),2)+pow((verrhopost+A),2))));
			//verseglength=sqrt(pow((zclass[jp]-zclass[ir]),2));
			///* verdelsquaredV in mV */
			//verdelsquaredV = (verVpre+verVpost-2*verVmid)*1000;
                        //if (verdelsquaredV > THRESHD2V)
			//    veraxonflagpos[nl]=1;
			//if (-verdelsquaredV > THRESHD2V)
			//    veraxonflagneg[nl]=1;

                        // Vertical axon segment  (IS Version)
                        veraxonflagpos[nl]=0;
			veraxonflagneg[nl]=0;
			verzmid=(2*zclass[ir]-2*ISLENGTH)/2;
                        verxmid=pre[0];
			verymid=pre[1];
			verrhomid=sqrt(pow((DISKX-verxmid),2)+pow((DISKY-verymid),2))*25;
			verVmid=(coeff)*asin(2*A/(sqrt(pow((STIMPOS-verzmid),2)+pow((verrhomid-A),2))+sqrt(pow((STIMPOS-verzmid),2)+pow((verrhomid+A),2))));
			verrhopre=sqrt(pow((DISKX-pre[0]),2)+pow((DISKY-pre[1]),2))*25;
			verVpre=(coeff)*asin(2*A/(sqrt(pow((STIMPOS-zclass[ir]),2)+pow((verrhopre-A),2))+sqrt(pow((STIMPOS-zclass[ir]),2)+pow((verrhopre+A),2))));
			verrhopost=sqrt(pow((DISKX-pre[0]),2)+pow((DISKY-pre[1]),2))*25;
			verVpost=(coeff)*asin(2*A/(sqrt(pow((STIMPOS-zclass[ir]+2*ISLENGTH),2)+pow((verrhopost-A),2))+sqrt(pow((STIMPOS-zclass[ir]+2*ISLENGTH),2)+pow((verrhopost+A),2))));
			verseglength=2*ISLENGTH; //sqrt(pow((zclass[jp]-zclass[ir]),2));
			/* verdelsquaredV in mV */
			verdelsquaredV = (verVpre+verVpost-2*verVmid)*1000;
                        if (verdelsquaredV > THRESHD2V)
			    veraxonflagpos[nl]=1;
			if (-verdelsquaredV > THRESHD2V)
			    veraxonflagneg[nl]=1;


			/* if (-verdelsquaredV > THRESHD2V){
			     currentinjsum+=(-verdelsquaredV)*Gi;
			     synapsesum=synapsesum+1;
			} */

			//if ( (horaxonflagpos[nl]==1) || (veraxonflagpos[nl]==1) )
			if ( veraxonflagpos[nl]==1 )
			    axonflag[nl]=1;

			//if ( (horaxonflagneg[nl]==1) || (veraxonflagneg[nl]==1) )
			if ( veraxonflagneg[nl]==1 )
			    axonflag[nl]=2;

                        ns+=zapisz2(nr,tab[nl],jp,ir,dist[nl],outputfile,graphfile,axflagfile,axonflag[nl]);
			/* printf("\npostindex=%d, postclass=%d, preindex=%d, preclass=%d",nr,ir+1,tab[nl],jp+1);*/  /* network linkage checks */
			// tab[nl]=-1;

			}
			
                      } /* koniec petli losowania */
                  } /* jp */

	}
     }
  fclose(outputfile);
  fclose(graphfile);
  fclose(axflagfile);

  // avecurrentinj=currentinjsum/synapsesum;
  //printf("currentinjsum=%f,avecurrentinj=%f,synapsesum=%i\n",currentinjsum,avecurrentinj,synapsesum);  /* network linkage checks */


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


//sort_delay  - uzywana w qsort() przy porzadkowaniu delay'ow

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 - uzywana w 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. Szczegoly generacji tablicy otoczenia:

*/

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)
{
 int i,j;
  do
    {
     i=randomm(index);
    }while(tablica[i]<0);

  j=tablica[i];
  tablica[i]=-1;
  return j;
}

int losuj2(int *tablica, double *distance, double *costht,double *sintht,int index, double Ap, double lambda, int classpost, int classpre,int *precoord)    
/* zwraca index w tab, a nie numer neuronu jak losuj() */
{
  int i,j,flagspatial;
  //double Ap=1.0;
  //double lambda=6;
  double los,losang,los1,los2,cutoff;
  double p3nearrad,p3farave,p3spotrad;
  double stripsize,preyy,postycoord,prefloor,postfloor;
  double precheck,postcheck;

  // i=randomm(index);

    do
    {
     i=randomm(index);
     // printf("i=%d,tablica[i]=%d,index=%d\n",i,tablica[i],index);  /* network linkage checks */
    }while(tablica[i]<0); 
  
  switch(classpre)
                     {

                        case 0:                           /* Pre - Layer II/III Pyramid */

			   /* Patch Handling */

			   p3nearrad=PIIINEARRAD/SPATIALSCALE;
                           stripsize=STRIPWIDTH/SPATIALSCALE;
			   //p3farave=PIIIFARAVE/SPATIALSCALE;
			   //p3spotrad=PIIISPOTRAD/SPATIALSCALE;
                           los=drand48();
                           // los2= Ap*exp(-1.0*distance[i]/lambda);
			   //losang=drand48();
                           // los2=1;
			   flagspatial=0;
			   if (distance[i]<=p3nearrad)
			       flagspatial=1;
                           else {
			       /* Stripflag A=even, B=odd */ 
                               preyy=precoord[1];
			       postycoord=posty[i];
                               prefloor=floor(preyy/stripsize);
                               postfloor=floor(postycoord/stripsize);
			       precheck=prefloor-2*floor(prefloor/2.0);
			       postcheck=postfloor-2*floor(postfloor/2.0);
                               //printf("precheck=%f,postcheck=%f,preyy=%f\n",precheck,postcheck,preyy); 
			       if ( precheck == postcheck)
			   	   flagspatial=0;
			       else
                                   flagspatial=0;

			       //printf("precheck=%f,postcheck=%f,flagspatial=%d\n",precheck,postcheck,flagspatial); 

                           }
			 
			   /* Angular isotropy routine */
			   //los1=(costht[i]*cos(PIIIANG)+sintht[i]*sin(PIIIANG)); 
			   //los1*=los1; 
			   //los1*=los1;
			   //los1*=los1; 
                           /*printf("p3nearrad=%lf,p3farave=%lf,p3spotrad=%lf,flagspatial=%i,distance=%lf\n",p3nearrad,p3farave,p3spotrad,flagspatial,distance[i]);*/  

                           if( ( flagspatial ) /*&&  (losang < los1)*/ ) 
                           {
                              return i;
                           }
                           else
                           {
			       // tablica[i]=-1; 
                              return -1;
                           }
                           break;

			   //los=drand48();
                           //// los2= Ap*exp(-1.0*distance[i]/lambda);
                           //los2=1;
                           ////printf("index=%i, neuron=%i, distance=%lf, drand48=%3.2lf, prob= %3.2lf\n",i,j,distance[i],los,los2);  
                           //if( los < los2 )
                           //{
                           //   return i;
                           //}
                           //else
                           //{
			   //    // tablica[i]=-1; 
                           //   return -1;
                           //}  
                           //break;


                        case 1:                           /* Pre - Layer IV Spiny */
                           los=drand48();
                           // los2= Ap*exp(-1.0*distance[i]/lambda);
                           los2=1;
                           //printf("index=%i, neuron=%i, distance=%lf, drand48=%3.2lf, prob= %3.2lf\n",i,j,distance[i],los,los2);  
                           if( los < los2 )
                           {
                              return i;
                           }
                           else
                           {
			       // tablica[i]=-1; 
                              return -1;
                           }  
                           break;
                        case 2:                           /* Pre - Layer V Pyramid */
                           los=drand48();
                           // los2= Ap*exp(-1.0*distance[i]/lambda);
			   //losang=drand48();
                           los2=1;
			   //los1=(costht[i]*cos(PVANG)+sintht[i]*sin(PVANG));
			   //los1*=los1;
			   //los1*=los1;
			   //los1*=los1;
                           //printf("index=%i, neuron=%i, distance=%lf, drand48=%3.2lf, prob= %3.2lf\n",i,j,distance[i],los,los2);  
			   if( ( los < los2 ) /*  && (losang < los1) */ )
                           {
                              return i;
                           }
                           else
                           {
			       // tablica[i]=-1; 
                              return -1;
                           }                           
                           break;
                        case 3:                           /* Pre - Layer VI Pyramid */
                           los=drand48();
                           // los2= Ap*exp(-1.0*distance[i]/lambda);
			   //losang=drand48();
                           los2=1;
			   //los1=(costht[i]*cos(PVIANG)+sintht[i]*sin(PVIANG));
			   //los1*=los1;
			   //los1*=los1;
			   //los1*=los1;
                           //printf("index=%i, neuron=%i, distance=%lf, drand48=%3.2lf, prob= %3.2lf\n",i,j,distance[i],los,los2);  
			   if( ( los < los2 ) /*  && (losang < los1) */ )
                           {
                              return i;
                           }
                           else
                           {
			       // tablica[i]=-1; 
                              return -1;
                           }                          
                           break;
                        case 4:                           /* Pre - Basket Cell */
			   cutoff = BASKETNEAR/SPATIALSCALE;
			   if ( distance[i] <= cutoff ){       /* Isotropic connections to 50 microns */
			       return i;
                           }
			   else { 
                              los=drand48();
			      //losang=drand48();
			      lambda = BASKETLAMBDA/SPATIALSCALE;
                              los2= Ap*exp(-1.0*distance[i]/lambda); 
			      //los1=(costht[i]*cos(BASKETANG)+sintht[i]*sin(BASKETANG));
			      //los1*=los1;
			      //los1*=los1;
			      //los1*=los1;
			      // los2=1;
                              /* printf("index=%i, neuron=%i, distance=%lf, drand48=%3.2lf, prob= %3.2lf\n",i,j,distance[i],los,los2); */ 
                              if( (los < los2) /* && (losang < los1) */ )
                              {
				  /* printf("index=%i, neuron=%i, costht=%lf, drand48=%3.2lf, prob= %3.2lf\n",i,j,costht[i],losang,los1);*/
                                 return i;
                              }
                              else
                              {
				  // tablica[i]=-1; 
                                 return -1;
                              }
                           }                           
                           break; 
                        case 5:                           /* Pre - Double Bouquet */
                           los=drand48();
                           // los2= Ap*exp(-1.0*distance[i]/lambda);
                           los2=1;
                           //printf("index=%i, neuron=%i, distance=%lf, drand48=%3.2lf, prob= %3.2lf\n",i,j,distance[i],los,los2);  
                           if( los < los2 )
                           {
                              return i;
                           }
                           else
                           {
			       // tablica[i]=-1; 
                              return -1;
                           }                           
                           break;  
                        case 6:                           /* Pre - Chandelier Cell  */
                           los=drand48();
                           // los2= Ap*exp(-1.0*distance[i]/lambda);
                           los2=1;
                           //printf("index=%i, neuron=%i, distance=%lf, drand48=%3.2lf, prob= %3.2lf\n",i,j,distance[i],los,los2);  
                           if( los < los2 )
                           {
                              return i;
                           }
                           else
                           {
			       // tablica[i]=-1; 
                              return -1;
                           }
                           break;
                        default:
                           fprintf(stderr,"Problem with neuron class in losuj2 %d\n",classpre);
                           exit(0);
                      }

}

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;
syn_tmp.neuron_pre=inr;
syn_tmp.weight=w[i][j]+rozrzut(sdw[i][j]);
syn_tmp.delay =d[i][j]+rozrzut(sdd[i][j]);
//printf("%i %i %i %i\n",syn_tmp.neuron_pre, syn_tmp.neuron_post , syn_tmp.weight,syn_tmp.delay );
return(fwrite(&syn_tmp,sizeof(struct SYN_TMP),1,outputfile));

}

int zapisz2(int inr, int inl, int i, int j, double distance,  FILE * outputfile,  FILE * graphfile, FILE * axflagfile, int axflag)
{
struct SYN_TMP syn_tmp;
extern  int rozrzut(int);
int axonal_del=5; /* was 5 */
syn_tmp.neuron_post=inl;
syn_tmp.neuron_pre=inr;
syn_tmp.preclass_num=j;
syn_tmp.weight=w[i][j]+rozrzut(sdw[i][j]);
syn_tmp.delay =d[i][j]+rozrzut(sdd[i][j]);
//printf("distance=%3.2lf, axonal del=%3.2lf, delay=%i, delay + axonal del=%i\n", distance, distance*axonal_del, syn_tmp.delay, syn_tmp.delay+(int)rint(distance*axonal_del));
syn_tmp.axdelay=(int)rint(distance*axonal_del)+3*rozrzut(sdd[i][j]);
// printf("%i %i %i %i %i\n",syn_tmp.neuron_pre, syn_tmp.neuron_post , syn_tmp.weight, syn_tmp.delay , syn_tmp.axdelay);
if ((j)==4){
 fwrite(&inr,sizeof(ushort),1,graphfile);
 fwrite(&j,sizeof(ushort),1,graphfile);
 fwrite(&inl,sizeof(ushort),1,graphfile);
 fwrite(&i,sizeof(ushort),1,graphfile);
}
/* printf("axflag= %i\n",axflag); */
 fwrite(&axflag,sizeof(ushort),1,axflagfile);

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[jp][ir];
   /* obliczamy najpierw wspolrzedne i[3] neuronu na siatce jp, ktory jest najblizej x,y,z */ /*???*/
  
  
  no=0;
  for(j=0;j<3;j++)
         i[j]=xyz[j]/k[jp][j];  

  /* teraz otoczenie
     wersja bez zawijania (brzegowe maja mniejsze otoczenie) */ /*???*/

  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){
              dist[no]=sqrt((ix-i[0])*(ix-i[0])+(iy-i[1])*(iy-i[1])+(iz-i[2])*(iz-i[2]));
              costheta[no]=(ix-i[0])/dist[no];
	      sintheta[no]=(iy-i[1])/dist[no];
	      /* printf("ix=%i, iy=%i, i[0]=%i,i[1]=%i,costheta=%lf,sintheta=%lf,dist=%lf\n",ix,iy,i[0],i[1],costheta[no],sintheta[no],dist[no]);*/
              // printf("%6.1lf,",dist[no]); /* tablica z odleglosciami - nowa wersja 01/29/03 */
	      postx[no]=ix;
	      posty[no]=iy;

              tab[no++]=j;}
          }
        }
     }
//printf("\n\n");
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;
   /* tab must be declared earlier with sufficient size */

   ot = o[jp][ir];
   /* obliczamy najpierw wspolrzedne i[3] neuronu na siatce jp, ktory jest najblizej x,y,z */ /* ??? */
  
  
  no=0;
  for(j=0;j<3;j++)
         i[j]=xyz[j]/k[jp][j];  

  /* Torus case */

      
 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;
}



Loading data, please wait...