TITLE Stochastic Hodgkin and Huxley model incorporating channel noise (approximate effective version).

COMMENT

This mod-file implementes a stochastic version of the HH model incorporating channel noise.
This version is the approximate ``effective'' version, i.e. it employs a Ornstein-Uhlenbeck process with
the same mean and variance of the Markov model. Note that the time constant of the exponentially
decaying covariance function of this effective model is just an approximation of the time constant
of the Markov model.

Author: Daniele Linaro - daniele.linaro@unige.it
Date: September 2010

ENDCOMMENT

UNITS {
    (mA) = (milliamp)
    (mV) = (millivolt)
    (S)  = (siemens)
    (pS) = (picosiemens)
    (um) = (micrometer)
} : end UNITS


NEURON {
    SUFFIX HHcnf
    USEION na READ ena WRITE ina
    USEION k READ ek WRITE ik
    NONSPECIFIC_CURRENT il
    RANGE el, gl
    RANGE gnabar, gkbar
    RANGE gna, gk
    RANGE gamma_na, gamma_k
    RANGE m_inf, h_inf, n_inf
    RANGE tau_m, tau_h, tau_n, tau_y, tau_z
    RANGE var_y, var_z
    RANGE noise_y, noise_z
    RANGE Nna, Nk
    RANGE seed    
    : these auxiliary variables are only needed if the exact method of solution is used    
    RANGE mu_y, mu_z
    THREADSAFE
} : end NEURON


PARAMETER {
    gnabar  = 0.12   (S/cm2)
    gkbar   = 0.036  (S/cm2)
    gl      = 0.0003 (S/cm2)  
    el       = -54.3 (mV)       : leakage reversal potential
    gamma_na = 20  (pS)		: single channel sodium conductance
    gamma_k  = 15  (pS)		: single channel potassium conductance
    seed = 5061983              : always use the same seed
} : end PARAMETER


STATE {
    m h n yy zz
} : end STATE


ASSIGNED {
    ina        (mA/cm2)
    ik         (mA/cm2)
    il         (mA/cm2)
    gna        (S/cm2)
    gk         (S/cm2)
    ena        (mV)
    ek         (mV)
    
    dt         (ms)
    area       (um2)
    celsius    (degC)
    v          (mV)
        
    Nna  (1) : number of potassium channels
    Nk   (1) : number of sodium channels
    
    m_inf h_inf n_inf
    noise_y noise_z
    var_y (ms2) var_z (ms2)
    tau_m (ms) tau_h (ms) tau_n (ms) tau_y (ms) tau_z (ms)
    : these auxiliary variables are only needed if the exact method of solution is used
    mu_y mu_z
} : end ASSIGNED

INITIAL {
    Nna = ceil(((1e-8)*area)*(gnabar)/((1e-12)*gamma_na))   : area in um2 -> 1e-8*area in cm2; gnabar in S/cm2; gamma_na in pS -> 1e-12*gamma_na in S
    Nk = ceil(((1e-8)*area)*(gkbar)/((1e-12)*gamma_k))   : area in um2 -> 1e-8*area in cm2; gkbar in S/cm2; gamma_k in pS -> 1e-12*gamma_k in S
    
    rates(v)
    m = m_inf
    h = h_inf
    n = n_inf
    yy = 0.
    zz = 0.
    set_seed(seed)
} : end INITIAL


BREAKPOINT {
    SOLVE states
    gna = gnabar * (m*m*m*h + zz)
    if (gna < 0) {
	gna = 0
    }
    gk = gkbar * (n*n*n*n + yy)
    if (gk < 0) {
	gk = 0
    }
    ina = gna * (v - ena)
    ik  = gk * (v - ek)
    il  = gl * (v - el)
} : end BREAKPOINT


PROCEDURE states() {
    rates(v)
    m = m + dt * (m_inf-m)/tau_m
    h = h + dt * (h_inf-h)/tau_h
    n = n + dt * (n_inf-n)/tau_n
    : Exact
    yy = yy*mu_y + noise_y
    zz = zz*mu_z + noise_z
    : Euler-Maruyama
    :yy = yy - dt * yy / tau_y + noise_y
    :zz = zz - dt * zz / tau_z + noise_z
    
    :consistency() 
    VERBATIM
    return 0;
    ENDVERBATIM
} : end PROCEDURE states()


PROCEDURE rates(Vm (mV)) { 
    LOCAL a,b,m3_inf,n4_inf,q10,sum,var_z1,var_z2,var_z3,var_z4,var_z5,var_z6,var_z7,tau_z1,tau_z2,tau_z3,tau_z4,tau_z5,tau_z6,tau_z7,var_y1,var_y2,var_y3,var_y4,one_minus_m,one_minus_h,one_minus_n
    
    UNITSOFF
    
    q10 = 3.^((celsius-6.3)/10.)
    
    : alpha_m and beta_m
    a = alpham(Vm)
    b = betam(Vm)
    sum = a+b
    tau_m = 1. / (q10*sum)
    m_inf = a / sum
    one_minus_m = 1. - m_inf
    m3_inf = m_inf*m_inf*m_inf

    : alpha_h and beta_h
    a = alphah(Vm)
    b = betah(Vm)
    sum = a+b
    tau_h = 1. / (q10*sum)
    h_inf = a / sum
    one_minus_h = 1. - h_inf
    
    tau_z1 = tau_h
    tau_z2 = tau_m
    tau_z3 = tau_m/2
    tau_z4 = tau_m/3
    tau_z5 = tau_m*tau_h/(tau_m+tau_h)
    tau_z6 = tau_m*tau_h/(tau_m+2*tau_h)
    tau_z7 = tau_m*tau_h/(tau_m+3*tau_h)
    var_z1 = 1.0 / Nna * m3_inf*m3_inf*h_inf * one_minus_h
    var_z2 = 3.0 / Nna * m3_inf*m_inf*m_inf*h_inf*h_inf * one_minus_m
    var_z3 = 3.0 / Nna * m3_inf*m_inf*h_inf*h_inf * one_minus_m*one_minus_m
    var_z4 = 1.0 / Nna * m3_inf*h_inf*h_inf * one_minus_m*one_minus_m*one_minus_m
    var_z5 = 3.0 / Nna * m3_inf*m_inf*m_inf*h_inf * one_minus_m*one_minus_h
    var_z6 = 3.0 / Nna * m3_inf*m_inf*h_inf * one_minus_m*one_minus_m*one_minus_h
    var_z7 = 1.0 / Nna * m3_inf*h_inf * one_minus_m*one_minus_m*one_minus_m*one_minus_h
    
    tau_z = (var_z1+var_z2+var_z3+var_z4+var_z5+var_z6+var_z7) / (var_z1/tau_z1+var_z2/tau_z2+var_z3/tau_z3+var_z4/tau_z4+var_z5/tau_z5+var_z6/tau_z6+var_z7/tau_z7)
    var_z = 1.0 / Nna * m3_inf * h_inf * (1.0 - m3_inf*h_inf)
     
    : Exact
    mu_z = exp(-dt/tau_z)
    noise_z = sqrt(var_z * (1-mu_z*mu_z)) * normrand(0,1)
    : Euler-Maruyama
    :noise_z = sqrt(2 * dt * var_z / tau_z) * normrand(0,1)
    
    : alpha_n and beta_n
    a = alphan(Vm)
    b = betan(Vm)
    sum = a+b
    tau_n = 1. / (q10*sum)
    n_inf = a / sum
    one_minus_n = 1. - n_inf
    n4_inf = n_inf * n_inf * n_inf * n_inf
    
    var_y1 = 4.0/Nk * n4_inf*n_inf*n_inf*n_inf * one_minus_n
    var_y2 = 6.0/Nk * n4_inf*n_inf*n_inf * one_minus_n*one_minus_n
    var_y3 = 4.0/Nk * n4_inf*n_inf * one_minus_n*one_minus_n*one_minus_n
    var_y4 = 1.0/Nk * n4_inf * one_minus_n*one_minus_n*one_minus_n*one_minus_n

    tau_y = (var_y1+var_y2+var_y3+var_y4) / (var_y1/tau_n + var_y2/tau_n/2 + var_y3/tau_n/3 + var_y4/tau_n/4)
    var_y = 1.0/Nk * n4_inf * (1.0-n4_inf)
    
    : Exact
    mu_y = exp(-dt/tau_y)
    noise_y = sqrt(var_y * (1-mu_y*mu_y)) * normrand(0,1)
    : Euler-Maruyama
    :noise_y = sqrt(2 * dt * var_y / tau_y) * normrand(0,1)
    
    UNITSON
}

PROCEDURE consistency(){
    if (m > 1) {
	m = 1
    }
    if (m < 0) {
	m = 0
    }
    if (h > 1) {
        h = 1
    }
    if (h < 0) {
	h = 0
    }
    if (n > 1) {
	n = 1
    }
    if (n < 0) {
	n = 0
    } 
}


FUNCTION vtrap(x,y) {  :Traps for 0 in denominator of rate eqns.
    if (fabs(x/y) < 1e-6) {
        vtrap = y*(1 - x/y/2)
    }else{
        vtrap = x/(exp(x/y) - 1)
    }
}


FUNCTION alpham(Vm (mV)) (/ms) {
    UNITSOFF
    alpham = .1 * vtrap(-(Vm+40),10)
    UNITSON
}


FUNCTION betam(Vm (mV)) (/ms) {
    UNITSOFF
    betam =  4 * exp(-(Vm+65)/18)
    UNITSON
}


FUNCTION alphah(Vm (mV)) (/ms) {
    UNITSOFF
    alphah = .07 * exp(-(Vm+65)/20)
    UNITSON
}


FUNCTION betah(Vm (mV)) (/ms) {
    UNITSOFF
    betah = 1 / (exp(-(Vm+35)/10) + 1)
    UNITSON
}


FUNCTION alphan(Vm (mV)) (/ms) {
    UNITSOFF
    alphan = .01*vtrap(-(Vm+55),10) 
    UNITSON
}


FUNCTION betan(Vm (mV)) (/ms) {
    UNITSOFF
    betan = .125*exp(-(Vm+65)/80)
    UNITSON
}

COMMENT
PROCEDURE evaluate_fct(v(mV)) { 
    LOCAL a, b, q10, sum
    UNITSOFF
    
    q10 = 3^((celsius-6.3)/10)
    
    : alpha_m and beta_m
    a = alpham(v)
    b = betam(v)
    sum = a+b
    tau_m = 1. / (q10*sum)
    m_inf = a / sum
    var_m = (a*(1-m) + b*m) / Nna
    m_noise = nom * sqrt(var_m * dt) * normrand(0,1)
        
    : alpha_h and beta_h
    a = alphah(v)
    b = betah(v)
    sum = a+b
    tau_h = 1. / (q10*sum)
    h_inf = a / sum
    var_h = (a*(1-h) + b*h) / Nna
    h_noise = noh * sqrt(var_h * dt) * normrand(0,1)
    
    : alpha_n and beta_n
    a = alphan(v)
    b = betan(v)
    sum = a+b
    tau_n = 1. / (q10*sum)
    n_inf = a / sum
    var_n = (a*(1-n) + b*n) / Nk
    n_noise = non * sqrt(var_n * dt) * normrand(0,1)
    UNITSON
}
ENDCOMMENT