Spiking neuron model of the basal ganglia (Humphries et al 2006)

 Download zip file 
Help downloading and running models
Accession:83559
A spiking neuron model of the basal ganglia (BG) circuit (striatum, STN, GP, SNr). Includes: parallel anatomical channels; tonic dopamine; dopamine receptors in striatum, STN, and GP; burst-firing in STN; GABAa, AMPA, and NMDA currents; effects of synaptic location. Model demonstrates selection and switching of input signals. Replicates experimental data on changes in slow-wave (<1 Hz) and gamma-band oscillations within BG nuclei following lesions and pharmacological manipulations.
Reference:
1 . Humphries MD, Stewart RD, Gurney KN (2006) A physiologically plausible model of action selection and oscillatory activity in the basal ganglia. J Neurosci 26:12921-42 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Realistic Network;
Brain Region(s)/Organism: Basal ganglia;
Cell Type(s): Neostriatum medium spiny direct pathway GABA cell; Subthalamus nucleus projection neuron; Globus pallidus neuron; Abstract integrate-and-fire leaky neuron;
Channel(s):
Gap Junctions:
Receptor(s): Dopaminergic Receptor;
Gene(s):
Transmitter(s): Dopamine; Gaba; Glutamate;
Simulation Environment: MATLAB;
Model Concept(s): Oscillations; Parkinson's; Action Selection/Decision Making; Sleep; Rebound firing;
Implementer(s): Humphries, Mark D [m.d.humphries at shef.ac.uk];
Search NeuronDB for information about:  Neostriatum medium spiny direct pathway GABA cell; Dopaminergic Receptor; Dopamine; Gaba; Glutamate;
/*=================================================================
 *
 * GHS_LIF_solver_shunt.C
 *
 *
 * Mark Humphries & Rob Stewart (6/1/2005)
 *=================================================================*/

#include <math.h>
#include "mex.h"

/* Input Arguments */
#define IN_N    prhs[0] 
#define IN_T    prhs[1]
#define LN_IN   prhs[2]
#define AMPA_IN   prhs[3]
#define NMDA_IN   prhs[4]
#define GABAA_IN   prhs[5]
#define LT_IN   prhs[6]
#define MAX_P   prhs[7]
#define MAX_S   prhs[8]
#define IG      prhs[9]
#define DELAYS  prhs[10]
#define BOUNDS  prhs[11]
#define N_PATHS prhs[12]
#define SPON    prhs[13]
#define T_ONSET prhs[14]
#define T_OFF   prhs[15]
#define P_CELLS prhs[16]
#define B_T1    prhs[17]
#define B_T2    prhs[18]
#define A_CA    prhs[19]
#define THETA_CA prhs[20]
#define C_CELLS prhs[21]
#define THETA   prhs[22]
#define MLIMIT  prhs[23]
#define TAU_AMPA   prhs[24]
#define TAU_NMDA   prhs[25]
#define TAU_GABAA   prhs[26]
#define TAU_M   prhs[27]
#define R_IN    prhs[28]

#define SCALAR_WS  prhs[29] /* weights = [SD1_w SD2_w STN_GPew STN_GPiw GPe_STNw GPe_GPiw GPe_GPew GPi_GPiw EXT_w STN_ext_ratio]; */
                            /* recall elements numbered from 0... For future expansion */

#define DPS     prhs[30]
#define IPS     prhs[31]
#define INIT_M  prhs[32]
#define INIT_I_GABAA  prhs[33]
#define INIT_I_SOMA  prhs[34]
#define INIT_I_PROX  prhs[35]
#define INIT_I_AMPA  prhs[36]
#define INIT_I_NMDA  prhs[37]

/* Output Arguments */
#define OUT_N       plhs[0]
#define OUT_T       plhs[1]
#define N_OUT       plhs[2]
#define N_EVENTS    plhs[3]
#define TRACE_VALS  plhs[4]
#define MEM_POT     plhs[5]
#define I_GABAA     plhs[6]
#define I_PROX      plhs[7]
#define I_SOMA      plhs[8]
#define I_AMPA      plhs[9]
#define I_NMDA      plhs[10]

/*-------------------- RANDOM NUMBER GENERATOR -----------------------------------*/
#define SHR3 (jz=jsr, jsr^=(jsr<<13), jsr^=(jsr>>17), jsr^=(jsr<<5),jz+jsr)
#define UNI (.5 + (signed) SHR3 * .2328306e-9)
#define RNOR (hz=SHR3, iz=hz&127, (abs(hz)<kn[iz])? hz*wn[iz] : nfix())
#define REXP (jz=SHR3, iz=jz&255, ( jz <ke[iz])? jz*we[iz] : efix())

static unsigned int iz,jz,jsr=123456789,kn[128],ke[256];
static int hz; static float wn[128],fn[128], we[256],fe[256];

float nfix(void) { /*provides RNOR if #define cannot */
    const float r = 3.442620f; static float x, y;
    for(;;){ x=hz*wn[iz];
        if(iz==0){ do{x=-log(UNI)*0.2904764; y=-log(UNI);} while(y+y<x*x);
        return (hz>0)? r+x : -r-x;
        }
        if( fn[iz]+UNI*(fn[iz-1]-fn[iz]) < exp(-.5*x*x) ) return x;
        hz=SHR3; iz=hz&127;if(abs(hz)<kn[iz]) return (hz*wn[iz]);
        } 
    }
    
float efix(void) { /*provides REXP if #define cannot */
    float x; 
    for(;;){
        if(iz==0) return (7.69711-log(UNI));
        x=jz*we[iz];
        if( fe[iz]+UNI*(fe[iz-1]-fe[iz]) < exp(-x) ) return (x);
        jz=SHR3; iz=(jz&255);
        if(jz<ke[iz]) return (jz*we[iz]);
    } 
}

/*--------This procedure sets the seed and creates the tables------*/
void zigset(unsigned int jsrseed) {
    const double m1 = 2147483648.0, m2 = 4294967296.;
    double dn=3.442619855899,tn=dn,vn=9.91256303526217e-3, q;
    double de=7.697117470131487, te=de, ve=3.949659822581572e-3;
    int i; jsr=jsrseed;

    /* Tables for RNOR: */ 
    q=vn/exp(-.5*dn*dn);
    kn[0]=(dn/q)*m1; kn[1]=0;
    wn[0]=q/m1; wn[127]=dn/m1;
    fn[0]=1.; fn[127]=exp(-.5*dn*dn);
    for(i=126;i>=1;i--) {
        dn=sqrt(-2.*log(vn/dn+exp(-.5*dn*dn)));
        kn[i+1]=(dn/tn)*m1; tn=dn;
        fn[i]=exp(-.5*dn*dn); wn[i]=dn/m1; 
    }
    
    /* Tables for REXP */ 
    q = ve/exp(-de);
    ke[0]=(de/q)*m2; ke[1]=0;
    we[0]=q/m2; we[255]=de/m2;
    fe[0]=1.; fe[255]=exp(-de);
    for(i=254;i>=1;i--) {
        de=-log(ve/de+exp(-de));
        ke[i+1]= (de/te)*m2; te=de;
        fe[i]=exp(-de); we[i]=de/m2; 
    }
}


/*------------------------------- SIMULATION CODE STARTS HERE -----------------------------------*/
static void build_ext( /* //doesn't look optimal - 04/02/04 */
                unsigned int    *extn,
                unsigned int    *extp,
                unsigned int    *n_ext,
                unsigned int    *n_paths,
                unsigned int    *delays,
                unsigned int    n_sources,
                unsigned int    p,
                unsigned int    tc,
                unsigned int    h,
                unsigned int    MaxEx
              )
{
    int path, newt;
    for(path = p; path < (n_paths[p]*n_sources+p); path+=n_sources){
      newt = (delays[path] + tc)%h; /* //circular addressing */
      extn[n_ext[newt]*h+newt] = p;      /* //store source */
      extp[n_ext[newt]*h+newt] = path;  /* //store path */
        n_ext[newt]++;
        if (n_ext[newt] >= MaxEx) mexErrMsgTxt("Bounding error. Increase MaxEx");
    }
}

static void IntFire3_GHS( 
                unsigned int    *out_n,
                unsigned int    *out_t,
                unsigned int    *n_out,
                unsigned int    *n_events,
                double          *trace_vals,

		double          *m,           /* initialises neuron state variables     */
		double          *gabaa,       /* and returns last values                */
		double          *ip,          /* so that network can be started          */
		double          *is,          /* using some snapshot                     */ 
		double          *ampa,        /* (excluding both refractory period and   */
		double          *nmda,        /* burst firing status)                    */

                unsigned int    *in_n,
                unsigned int    *in_t,
                mxArray         *LN,
                mxArray         *L_AMPA,
                mxArray         *L_NMDA,
                mxArray         *L_GABAA,
                mxArray         *LT,
                double          *max_prox,
                double          *max_soma,    
                double          *i_gate,
                unsigned int    *delays,
                unsigned int    *bounds,
                unsigned int    *n_paths,
                double          *spon,
                unsigned int    *t_onset,
                unsigned int    *t_off,
                unsigned int    *p_cells,
                unsigned int    *burst_t1, 
                unsigned int    *burst_t2, 
                double          *alphaCA, 
                double          *thetaCA,
                unsigned int    *ca_cells,
                double          *theta,
                double          *mlimit,
                double          *tau_AMPA,
                double          *tau_NMDA,
                double          *tau_GABAa,
                double          *tau_m,
                double          *R,
                double          *Wgts,
                double          *dps,
                unsigned int    *ips,
		double          *init_mem_pot, 
		double          *init_i_gabaa, 
		double          *init_i_soma, 
		double          *init_i_prox, 
		double          *init_i_ampa, 
		double          *init_i_nmda
              )
{
    /* unbundle parameters */
    double sigma = dps[0];     /* //Noise std */
    double ref = dps[1];       /* //Refractory membrane potential reset value */
    double dt = dps[2];
    double step_size = dps[3];
    double PSP_sigma = dps[4];
    double M_hat = dps[5];
    
    unsigned int h = ips[0];
    unsigned int MaxEx = ips[1];
    unsigned int MaxOut = ips[3];
    unsigned int n_neurons = ips[4];
    unsigned int n_sources = ips[5];
    unsigned int n_in = ips[6];
    unsigned int time_steps = ips[7];
    unsigned int ref_period = ips[8];
    unsigned int trace = ips[9];
    
    unsigned int n_ext_total = 0;
    unsigned int n_out_val = 0;
    unsigned int tc = 0;  
    unsigned int in = 0;
    unsigned int t,p,loop,loop2,source,first,last;
    int idum[] = {-1};
    double w, psp_noise, pre_m, soma_gate, prox_gate, total_gate;
    
    unsigned int p_index = 0;
    
    /* mxArray *AMPA, *NMDA, *GABAA, *IP, *IS, *M, */
    mxArray  *I_EXT, *I_CA, *B_END, *T_CA, *T1, *E_AMPA, *E_NMDA, *E_GABAA, *EM, *AS, *EXTN, *EXTP, *N_EXT;
    /* double *ampa, *nmda, *gabaa, *m, *ip, *is; */
    double *i_ext, *i_ca, *as, *e_ampa, *e_nmda, *e_gabaa, *em, *lw_ampa, *lw_nmda, *lw_gabaa;
    unsigned int *extn, *extp, *n_ext, *ln, *lt;
    unsigned int type;
    int *t1, *t_ca, *b_end;
    
    int t1_dims[] = {1,0};
    int en_dims[] = {1,0};
    int ne_dims[] = {1,0};

    t1_dims[1] = n_neurons;
    en_dims[0] = h; en_dims[1] = MaxEx;
    ne_dims[1] = h;
    
    /*Create state variable arrays*/
    /* AMPA = mxCreateDoubleMatrix(1,n_neurons,mxREAL);
    NMDA = mxCreateDoubleMatrix(1,n_neurons,mxREAL);
    GABAA = mxCreateDoubleMatrix(1,n_neurons,mxREAL);
    IP = mxCreateDoubleMatrix(1,n_neurons,mxREAL);
    IS = mxCreateDoubleMatrix(1,n_neurons,mxREAL);
    M = mxCreateDoubleMatrix(1,n_neurons,mxREAL); */

    I_EXT = mxCreateDoubleMatrix(1,n_neurons,mxREAL);
    I_CA = mxCreateDoubleMatrix(1,n_neurons,mxREAL);
    AS = mxCreateDoubleMatrix(1,n_neurons,mxREAL);
    E_AMPA = mxCreateDoubleMatrix(1,n_neurons,mxREAL);
    E_NMDA = mxCreateDoubleMatrix(1,n_neurons,mxREAL);
    E_GABAA = mxCreateDoubleMatrix(1,n_neurons,mxREAL);
    EM = mxCreateDoubleMatrix(1,n_neurons,mxREAL);
    B_END = mxCreateNumericArray(2,t1_dims,mxINT32_CLASS,mxREAL);
    T_CA = mxCreateNumericArray(2,t1_dims,mxINT32_CLASS,mxREAL);
    T1 = mxCreateNumericArray(2,t1_dims,mxINT32_CLASS,mxREAL);
    
    /*assign pointers to state variables*/
    /* ampa = mxGetPr(AMPA);
    nmda = mxGetPr(NMDA);
    gabaa = mxGetPr(GABAA);
    ip = mxGetPr(IP);
    is = mxGetPr(IS);
    m = mxGetPr(M); */
    i_ext = mxGetPr(I_EXT);
    i_ca = mxGetPr(I_CA);
    as = mxGetPr(AS);
    e_ampa = mxGetPr(E_AMPA);
    e_nmda = mxGetPr(E_NMDA);
    e_gabaa = mxGetPr(E_GABAA);
    em = mxGetPr(EM);
    b_end = (int*)mxGetData(B_END);
    t_ca = (int*)mxGetData(T_CA);
    t1 = (int*)mxGetData(T1);
    
    /*Create queue variables*/
    EXTN = mxCreateNumericArray(2, en_dims, mxUINT32_CLASS, mxREAL);
    EXTP = mxCreateNumericArray(2, en_dims, mxUINT32_CLASS, mxREAL);
    N_EXT = mxCreateNumericArray(2, ne_dims, mxUINT32_CLASS, mxREAL);
    
    /*assign pointers to queue variables*/
    extn = mxGetData(EXTN);  
    extp = mxGetData(EXTP);
    n_ext = mxGetData(N_EXT);
    
     /*----- set random number seed ---------*/
    zigset(1);
    
    /* //Initialise t1, ee, ei, em, as */
    for(p = 0; p < n_neurons; p++){
        t_ca[p] = -1;
        t1[p] = -1;
        e_ampa[p] = exp(-dt/tau_AMPA[p]);
        e_nmda[p] = exp(-dt/tau_NMDA[p]);
        e_gabaa[p] = exp(-dt/tau_GABAa[p]);
        em[p] = exp(-dt/tau_m[p]);
        as[p] = R[p]*(1-em[p]);
        i_ca[p] = 0;
        i_ext[p] = 0;
        b_end[p] = burst_t1[p] + burst_t2[p]; /* // sum these here for efficiency */
        /* // initialise membrane and synaptic variables (history) */
        m[p] = init_mem_pot[p]; /* RNOR * 0.002 + 0.00; */
        ip[p] = init_i_prox[p];
        is[p] = init_i_soma[p];
        gabaa[p] = init_i_gabaa[p]; /* RNOR * 0.0001 - 0.0001; */
        ampa[p] = init_i_ampa[p];  /* RNOR * 0.0001 + 0.0001; */
        nmda[p] = init_i_nmda[p];  /* RNOR * 0.0001 + 0.0001; */
        /* //printf("Initial is %f\n",m[p]); */	
     }    
    
    
    /*Begin main loop over time steps*/
    
    printf("Tracing neuron... %d \n",trace);
    for(t = 0; t < time_steps; t++){
        /* //first trace update 'trace' is theneuron index you specified in matlab */
        /* // *************** customise detailed traces here **************** */
        trace_vals[t] = gabaa[trace]; /* distal inhibtory */
        trace_vals[time_steps + t] = ampa[trace] + nmda[trace]; /* total excitatory current */
        trace_vals[time_steps*2 + t] = m[trace];                /* membrane potential */
        /* trace_vals[time_steps + t] = gabaa[trace]; ** distal inhibitory */
        /* trace_vals[time_steps + t] = ip[trace]; ** proximal inhibitory */
        /* trace_vals[time_steps + t] = is[trace]; ** somatic inhibitory */
        /* trace_vals[time_steps + t] = i_ca[trace]; ** Ca++ inhibitory */
        
      for(p = 0; p < n_neurons; p++){ /* //loop over all neurons */
	
	    pre_m = m[p];	/* // membrane potential on previous time-step was */
        	
            if ((t-t1[p]) > ref_period) { 
	      /* //printf("Somatic gate %f \n",1 - (is[p] / max_soma[p])); */
                
                /* Option 1: scale gate size by max for that neuron */
                /* soma_gate = 1 - (is[p] / max_soma[p]);
                prox_gate = 1 - (ip[p] / max_prox[p]);  */
                
                /* Option 2: scale gate size by max for the network */
                soma_gate = 1 - (is[p] / M_hat);
                prox_gate = 1 - (ip[p] / M_hat);
                total_gate = 1 - (prox_gate + soma_gate) / 2; 

                /* // hard-limit gates */
                if (soma_gate < 0) soma_gate = 0;
                if (prox_gate < 0) prox_gate = 0;
                if (total_gate < 0) total_gate = 0;
                
                /* trace gate value */
        if (p == trace) trace_vals[time_steps + t] = total_gate;  /* // trace gating variables */
                
                
		/* // ********** The membrane equation ****************** */
        
                /* // gate inputs and spontaneous current */
                /* //m[p] = m[p]*em[p] + ((((e[p] + i[p]) * prox_gate)  + spon[p] + i_ca[p]) * soma_gate) + i_ext[p]; // gate all psps */
                
                /* // no gating, just increase size of clamping current  */
                /* m[p] = m[p]*em[p] + ampa[p] + nmda[p] + gabaa[p] + spon[p] + i_ca[p] + i_ext[p] + i_gate[p] * total_gate;
                
                /* gating and clamp current combined - shunt I_ca */
               /* m[p] = m[p]*em[p] + (((ampa[p] + nmda[p] + gabaa[p]) * prox_gate)  + i_ca[p]) * soma_gate + i_ext[p] + spon[p] + i_gate[p] * total_gate; */

	       /* gating and clamp current combined - don't shunt shunt I_ca */
               m[p] = m[p]*em[p] + (((ampa[p] + nmda[p] + gabaa[p]) * prox_gate)) * soma_gate + i_ext[p] + spon[p] + i_ca[p] + i_gate[p] * total_gate;
                
                /* // add noise - do this here because now have two thresholds etc */
                m[p] = m[p] + RNOR * sigma;
                
                /* // hard-limit to replicate reversal potential  */
                if(m[p] < mlimit[p]) m[p] = mlimit[p];
            }
            
            /* // check for external current input */
            if (p_cells[p]==1 & t > t_onset[p_index] & t < t_off[p_index]) {
	      /* // printf("p_cell %f\n",p_cells[p]); */
                i_ext[p] = step_size * as[p];
            }
            else i_ext[p] = 0;
            
            /* // do Ca current checks - is burst cell, not bursting, and rebounds */
            if (ca_cells[p]==1 & i_ca[p]==0 & pre_m < thetaCA[p] & m[p] > thetaCA[p]){    /* // rebound bursting */
	      /* //printf("burst triggered at %d by pre m %f and m %f \n",t,pre_m,m[p]); */
                  t_ca[p] = 1;
                  i_ca[p] = alphaCA[p];
            }
            else{
                if (ca_cells[p]==1 & i_ca[p] > 0){
                      t_ca[p]++;
                      /* //printf("t CA %d\n",t_ca[p]); */
                      if (t_ca[p] >= burst_t1[p] & t_ca[p] <= b_end[p]){
			/* // calculation for ramp - move constant calculation to preamble for efficiency */
                          i_ca[p] = alphaCA[p] - (alphaCA[p] / burst_t2[p])*(t_ca[p] - burst_t1[p]);
                      }
                      else{
                          if(t_ca[p] > b_end[p]){  
                             t_ca[p] = 0;
                             i_ca[p] = 0;
                          }
                      }
                }
            }        
            /* // update PSP totals */
            gabaa[p] *= e_gabaa[p]; 
            ip[p] *= e_gabaa[p];
            is[p] *= e_gabaa[p];
            ampa[p] *= e_ampa[p];
            nmda[p] *= e_nmda[p];

            
            /* // ******************** end membrane potential ****************  */
        }
        
      
            
      /* //Process external events */
        for(loop = tc; loop < (n_ext[tc]*h+tc); loop+=h){
            source = extn[loop];
            ln = mxGetData(mxGetCell(LN,source));   /* //get pointer to target neurons */
            lw_ampa = mxGetPr(mxGetCell(L_AMPA,source));     /* //get pointer to connection weights */            
            lw_nmda = mxGetPr(mxGetCell(L_NMDA,source));     /* //get pointer to connection weights */            
            lw_gabaa = mxGetPr(mxGetCell(L_GABAA,source));     /* //get pointer to connection weights */            
            lt = mxGetData(mxGetCell(LT,source));     /* //get pointer to connection type */
            first = bounds[extp[loop]];
            last = bounds[extp[loop]+n_sources];
            for(loop2 = first; loop2 < last; loop2++){
                n_ext_total++;
                p = ln[loop2];      /* //Neuron */
                if (p != 0) {       /* // then is connnected */
			        type = lt[loop2];   /* //Type                 */
			        psp_noise =  RNOR * PSP_sigma * w;	/* // scale PSP noise by weight */
			        /* //psp_noise = 0; */
              
					  switch(type){
					    case 0: ampa[p] += lw_ampa[loop2] + psp_noise;          /* //excitatory event */
					            nmda[p] += lw_nmda[loop2] + psp_noise;
                                            break;
					    case 1: gabaa[p] += lw_gabaa[loop2] + psp_noise;          /* //inhibitory event - distal */
                                            break;
					    case 2: ip[p] += lw_gabaa[loop2] + psp_noise;         /* //inhibitory event - proximal */
                                            break;
					    case 3: is[p] += lw_gabaa[loop2] + psp_noise;         /* //inhibitory event - soma */
                                            break;
                        default: printf("Shouldn't be here \n");
                      }
                }                
            } /* //loop2     */
        } /* //loop - external events */
        
        /* // check if current has been injected - if so, increment time array */
        if (t>0 & fmod(t,t_off[p_index]) == 0){
	  /* //printf("On loop %d t_onset is %f, t_off %f\n",t,t_onset[p_index],t_off[p_index]); */
            p_index++;
        }
        
        n_ext[tc] = 0; /* //reset number of external events */
        /* //printf("PSP noise.. %f \n",psp_noise); */
       
        for(p = 0; p < n_neurons; p++){ /* //fire loop */
	  /* //if (m[p] > (theta[p] + (RNOR * sigma ))){ //fire if above threshold + noise */
	  if (m[p] > theta[p]) { /* //fire if above threshold - noise now added to membrane potential */
                if (p == trace)
		  trace_vals[time_steps*2 + t] = m[trace]+0.04; /* //fake spikes */
                m[p] = ref;
                t1[p] = t;
                out_n[n_out_val] = p+1;
                out_t[n_out_val] = t+1;
                n_out_val++;
                if (n_out_val > MaxOut) mexErrMsgTxt("Bounding error. Increase MaxOut");
                build_ext(extn,extp,n_ext,n_paths,delays,n_sources,p,tc,h,MaxEx);
            }   
        } /* //fire loop */

        /* //Add extrinsic events to queue */
        while ((in_t[in]) == t){
            build_ext(extn,extp,n_ext,n_paths,delays,n_sources,in_n[in],tc,h,MaxEx);
            if (++in >= n_in) break;
        }
        
        tc++;
        if (tc == h) tc = 0; /* //circular addressing */
    }
    n_out[0] = n_out_val;
    n_events[0] = n_ext_total;
    printf("number of spike sources: \t%d\n",n_sources);
    printf("number of events: \t\t\t%d\n",n_ext_total);
    printf("number of spikes: \t\t\t%d\n",n_out_val);
}

void mexFunction( int nlhs, mxArray *plhs[], 
          int nrhs, const mxArray *prhs[] )
{
    double *dps, *spon, *alphaCA, *thetaCA, *theta, *mlimit, *tau_AMPA, *tau_NMDA, *tau_GABAa, *tau_m, *R, *trace_vals, *max_prox, *max_soma, *i_gate, *Wgts, *mem_pot, *i_gabaa, *i_soma, *i_prox, *i_ampa, *i_nmda, *init_mem_pot, *init_i_gabaa, *init_i_soma, *init_i_prox, *init_i_ampa, *init_i_nmda;
    unsigned int *ips, *in_n, *in_t, *out_n, *out_t, *n_out, *n_events, *delays, *bounds, *n_paths,*p_cells,*t_onset, *t_off, *burst_t1, *burst_t2, *ca_cells;
    int dims[] = {1,0};
    int dims2[] = {1,1};
    
    /* Check for proper number of arguments */
    /*if (nrhs != 31) { 
        mexErrMsgTxt("Thirty-one input arguments required."); 
    } else if (nlhs > 5) {
        mexErrMsgTxt("Too many output arguments."); 
    } */

    /* Assign input pointers - pass cell arrays unprocessed */
    in_n = mxGetData(IN_N);
    in_t = mxGetData(IN_T);
    max_prox = mxGetPr(MAX_P);
    max_soma = mxGetPr(MAX_S);
    i_gate = mxGetPr(IG);
    delays = mxGetData(DELAYS);
    bounds = mxGetData(BOUNDS);
    n_paths = mxGetData(N_PATHS);
    spon = mxGetPr(SPON);
    t_onset = mxGetData(T_ONSET);
    t_off = mxGetData(T_OFF);
    p_cells = mxGetData(P_CELLS);
    burst_t1 = mxGetData(B_T1);
    burst_t2 = mxGetData(B_T2);
    alphaCA = mxGetPr(A_CA);
    thetaCA = mxGetPr(THETA_CA);
    ca_cells = mxGetData(C_CELLS);
    theta = mxGetPr(THETA);
    mlimit = mxGetPr(MLIMIT);
    tau_AMPA = mxGetPr(TAU_AMPA);
    tau_NMDA = mxGetPr(TAU_NMDA);
    tau_GABAa = mxGetPr(TAU_GABAA);
    tau_m = mxGetPr(TAU_M);
    R = mxGetPr(R_IN);
    Wgts = mxGetPr(SCALAR_WS);
    dps = mxGetPr(DPS);
    ips = mxGetData(IPS);
    init_mem_pot = mxGetData(INIT_M);
    init_i_gabaa = mxGetData(INIT_I_GABAA);
    init_i_soma = mxGetData(INIT_I_SOMA); 
    init_i_prox = mxGetData(INIT_I_PROX);  
    init_i_ampa = mxGetData(INIT_I_AMPA); 
    init_i_nmda = mxGetData(INIT_I_NMDA);


    /* Create matrices for the return argument */ 
    dims[1] = ips[3]; /* //MaxOut */
    OUT_N = mxCreateNumericArray(2, dims, mxUINT32_CLASS, mxREAL);
    OUT_T = mxCreateNumericArray(2, dims, mxUINT32_CLASS, mxREAL);
    N_OUT = mxCreateNumericArray(2, dims2, mxUINT32_CLASS, mxREAL);
    N_EVENTS = mxCreateNumericArray(2, dims2, mxUINT32_CLASS, mxREAL);
    TRACE_VALS = mxCreateDoubleMatrix(ips[7],3,mxREAL);
    MEM_POT = mxCreateDoubleMatrix(ips[4],1,mxREAL);    
    I_GABAA =  mxCreateDoubleMatrix(ips[4],1,mxREAL); 
    I_PROX =  mxCreateDoubleMatrix(ips[4],1,mxREAL); 
    I_SOMA =  mxCreateDoubleMatrix(ips[4],1,mxREAL); 
    I_AMPA =  mxCreateDoubleMatrix(ips[4],1,mxREAL); 
    I_NMDA =  mxCreateDoubleMatrix(ips[4],1,mxREAL); 

    /* Assign output pointers */ 
    out_n = mxGetData(OUT_N);
    out_t = mxGetData(OUT_T);
    n_out = mxGetData(N_OUT);
    n_events = mxGetData(N_EVENTS);
    trace_vals = mxGetPr(TRACE_VALS);
    mem_pot = mxGetData(MEM_POT);
    i_gabaa = mxGetData(I_GABAA);
    i_prox = mxGetData(I_PROX);
    i_soma = mxGetData(I_SOMA);
    i_ampa = mxGetData(I_AMPA);    
    i_nmda = mxGetData(I_NMDA);    


    

    /* Do the actual computations in a subroutine */
    IntFire3_GHS(out_n,out_t,n_out,n_events,trace_vals,mem_pot,i_gabaa,i_soma,i_prox,i_ampa,i_nmda,
		 in_n,in_t,LN_IN,AMPA_IN,NMDA_IN,GABAA_IN,LT_IN,max_prox,max_soma,i_gate,delays,
		 bounds,n_paths,spon,t_onset,t_off,p_cells,burst_t1,burst_t2,alphaCA,thetaCA,
		 ca_cells,theta,mlimit,tau_AMPA,tau_NMDA,tau_GABAa,tau_m,R,Wgts,dps,ips,
		 init_mem_pot, init_i_gabaa, init_i_prox, init_i_soma, init_i_ampa, init_i_nmda); 
}

Loading data, please wait...