theory_EIFfpt.c

/* 
[f0] = theory_EIFfpt(params,x0,mu_in,xi,r0,P0,freq)
Calculates the first passage time density of EIF, which can be used to calculate the single-neuron spike train power spectrum.
This cannot be done for models with voltage-activated conductances that make the spike train non-renewal.
The SDE is:
C*V' = gL(VL-V)+gL*Delta*exp((V-VT)/Delta)+mu+gL*sigma*sqrt(2*C/gL)*gaussianwhitenoise(t);
with reset at Vr, hard threshold at Vth, lower boundary at Vlb and linear decay of the membrane potential from threshold to reset

freq is a vector of frequencies 
 */

#include "mex.h"
#include "math.h"
#include "complex.h"

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{

int f, k, Nloop, Nfreq, kre, m;
double sigma2, mu_in, gL, C, Delta, VT, VL, Vth, Vlb, dV, Vr, tref, V, w, pi;
double *params, *freq, *f0r, *f0i;
double df, fmax;

/* get inputs */
params = mxGetPr(prhs[0]);
mu_in = mxGetScalar(prhs[1]);
sigma2 = mxGetScalar(prhs[2]);
freq = mxGetPr(prhs[3]);
Nfreq = mxGetN(prhs[3]);
m = mxGetM(prhs[3]);

if (Nfreq==1 && m!=1) {
    Nfreq = m;
    m = 1;
}

// df = mxGetPr(prhs[3]);
// fmax = mxGetPr(prhs[4]);
// 
// Nfreq = (int)(*fmax)/(*df)*2;
// freq = mxMalloc(Nfreq*sizeof(double));
// freq[0] = -(*fmax);
// for (f=1; f<Nfreq; f++)
//     freq[f] += (*df);

gL = params[0];
C = params[1];
Delta = params[2];
VT = params[3];
VL = params[4];
Vth = params[5];
Vlb = params[6];
dV = params[7];
Vr = params[8];
tref = params[9];

Nloop = (int)floor((Vth-Vlb)/dV); /* number of bins */
kre = (int)round((Vr-Vlb)/dV); /* index of reset potential */
pi=3.14159265;

/* allocate outputs */
plhs[0] = mxCreateDoubleMatrix(Nfreq,1,mxCOMPLEX);

f0r = mxGetPr(plhs[0]);
f0i = mxGetPi(plhs[0]);

double *I0, *G, *alpha, *beta;
double _Complex *pf, *po, *jf, *jo;

pf=mxMalloc(Nloop*sizeof(double _Complex));
po=mxMalloc(Nloop*sizeof(double _Complex));
jf=mxMalloc(Nloop*sizeof(double _Complex));
jo=mxMalloc(Nloop*sizeof(double _Complex));
I0=mxMalloc(Nloop*sizeof(double));
G=mxMalloc(Nloop*sizeof(double));
alpha=mxMalloc(Nloop*sizeof(double));
beta=mxMalloc(Nloop*sizeof(double));

/* initialize arrays for integration */

V=Vlb;
for (k=0; k<Nloop; k++){
     
    I0[k] = gL*(VL-V)+mu_in+gL*Delta*exp((V-VT)/Delta);
    G[k] = -I0[k]/(gL*sigma2);
    alpha[k] = exp(dV*G[k]);
    if (G[k]==0) beta[k] = 1/sigma2;
    else beta[k] = (alpha[k]-1)/(G[k]*sigma2);
    
    jf[k] = 0;
    jo[k] = 0;
    pf[k] = 0;
    po[k] = 0;
    
    V += dV;
}


for (f=0; f<Nfreq; f++) {
    
/* initialize flux and probability to 0 */ 
jf[Nloop-1] = 1+0*I; // boundary condition
jo[Nloop-1] = 0+0*I;
pf[Nloop-1] = 0+0*I;
po[Nloop-1] = 0+0*I;
    
w = freq[f]*2*pi;

V = Vth;

for (k = Nloop-2; k>=0; k--){
       
    /* boundary conditions */
    jf[k] = jf[k+1]+dV*I*w*pf[k+1];
    pf[k] = pf[k+1]*alpha[k+1]+C/gL*jf[k+1]*beta[k+1];
    
    /* inhomogenous component - current modulations */
    jo[k] = jo[k+1]+dV*I*w*po[k+1]; if (k==kre) jo[k] -= cexp(-I*w*tref);
    po[k] = po[k+1]*alpha[k+1]+C/gL*jo[k+1]*beta[k+1];
    
    V -= dV;
    
}

f0r[f] = creal(-jo[0]/jf[0]);
f0i[f] = cimag(-jo[0]/jf[0]);

} // frequency loop

mxFree(pf);
mxFree(po);
mxFree(jf);
mxFree(jo);
mxFree(I0);
mxFree(G);
mxFree(alpha);
mxFree(beta);

} //mexFunction