Accurate and fast simulation of channel noise in conductance-based model neurons (Linaro et al 2011)

 Download zip file   Auto-launch 
Help downloading and running models
Accession:127992
We introduce and operatively present a general method to simulate channel noise in conductance-based model neurons, with modest computational overheads. Our approach may be considered as an accurate generalization of previous proposal methods, to the case of voltage-, ion-, and ligand-gated channels with arbitrary complexity. We focus on the discrete Markov process descriptions, routinely employed in experimental identification of voltage-gated channels and synaptic receptors.
Reference:
1 . Linaro D, Storace M, Giugliano M (2011) Accurate and fast simulation of channel noise in conductance-based model neurons by diffusion approximation. PLoS Comput Biol 7:e1001102 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Neuron or other electrically excitable cell;
Brain Region(s)/Organism: Neocortex;
Cell Type(s): Neocortex U1 L2/6 pyramidal intratelencephalic GLU cell; Neocortex U1 L5B pyramidal pyramidal tract GLU cell;
Channel(s): I Na,t; I K;
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: NEURON; C or C++ program; Python;
Model Concept(s): Ion Channel Kinetics; Simplified Models; Methods; Markov-type model;
Implementer(s): Linaro, Daniele [daniele.linaro at unige.it];
Search NeuronDB for information about:  Neocortex U1 L5B pyramidal pyramidal tract GLU cell; Neocortex U1 L2/6 pyramidal intratelencephalic GLU cell; I Na,t; I K;
/
HHcn
mod-files
HHcn.mod
HHcnf.mod
HHmicro.mod
                            
TITLE Stochastic Hodgkin and Huxley model incorporating channel noise (effective version).

COMMENT

This mod-file implementes a stochastic version of the HH model incorporating channel noise.
This version is the ``effective'' version, i.e. it employs a Ornstein-Uhlenbeck process with
the same stochastic properties (mean and covariance) as 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 HHcn
    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
    RANGE tau_y1, tau_y2, tau_y3, tau_y4
    RANGE tau_z1, tau_z2, tau_z3, tau_z4, tau_z5, tau_z6, tau_z7
    RANGE var_y1, var_y2, var_y3, var_y4
    RANGE var_z1, var_z2, var_z3, var_z4, var_z5, var_z6, var_z7
    RANGE noise_y1, noise_y2, noise_y3, noise_y4
    RANGE noise_z1, noise_z2, noise_z3, noise_z4, noise_z5, noise_z6, noise_z7
    RANGE Nna, Nk
    RANGE seed    
    : these auxiliary variables are only needed if the exact method of solution is used    
    RANGE mu_y1, mu_y2, mu_y3, mu_y4
    RANGE mu_z1, mu_z2, mu_z3, mu_z4, mu_z5, mu_z6, mu_z7
    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 = 10  (pS)		: single channel sodium conductance
    gamma_k  = 10  (pS)		: single channel potassium conductance
    seed = 5061983              : always use the same seed
} : end PARAMETER


STATE {
    m h n
    y1 y2 y3 y4
    z1 z2 z3 z4 z5 z6 z7
} : 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_y1 noise_y2 noise_y3 noise_y4
    noise_z1 noise_z2 noise_z3 noise_z4 noise_z5 noise_z6 noise_z7
    var_y1 (ms2) var_y2 (ms2) var_y3 (ms2) var_y4 (ms2) 
    var_z1 (ms2) var_z2 (ms2) var_z3 (ms2) var_z4 (ms2) var_z5 (ms2) var_z6 (ms2) var_z7 (ms2)
    tau_m (ms) tau_h (ms) tau_n (ms)
    tau_y1 (ms) tau_y2 (ms) tau_y3 (ms) tau_y4 (ms) 
    tau_z1 (ms) tau_z2 (ms) tau_z3 (ms) tau_z4 (ms) tau_z5 (ms) tau_z6 (ms) tau_z7 (ms)
    : these auxiliary variables are only needed if the exact method of solution is used
    mu_y1 mu_y2 mu_y3 mu_y4
    mu_z1 mu_z2 mu_z3 mu_z4 mu_z5 mu_z6 mu_z7
} : 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
    y1 = 0.
    y2 = 0.
    y3 = 0.
    y4 = 0.
    z1 = 0.
    z2 = 0.
    z3 = 0.
    z4 = 0.
    z5 = 0.
    z6 = 0.
    z7 = 0.
    set_seed(seed)
} : end INITIAL


BREAKPOINT {
    SOLVE states
    gna = gnabar * (m*m*m*h + z1+z2+z3+z4+z5+z6+z7)
    if (gna < 0) {
	gna = 0
    }
    gk = gkbar * (n*n*n*n + y1+y2+y3+y4)
    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
    y1 = y1*mu_y1 + noise_y1
    y2 = y2*mu_y2 + noise_y2
    y3 = y3*mu_y3 + noise_y3
    y4 = y4*mu_y4 + noise_y4
    z1 = z1*mu_z1 + noise_z1
    z2 = z2*mu_z2 + noise_z2
    z3 = z3*mu_z3 + noise_z3
    z4 = z4*mu_z4 + noise_z4
    z5 = z5*mu_z5 + noise_z5
    z6 = z6*mu_z6 + noise_z6
    z7 = z7*mu_z7 + noise_z7
    : Euler-Maruyama
    :y1 = y1 - dt * y1 / tau_y1 + noise_y1
    :y2 = y2 - dt * y2 / tau_y2 + noise_y2
    :y3 = y3 - dt * y3 / tau_y3 + noise_y3
    :y4 = y4 - dt * y4 / tau_y4 + noise_y4
    :z1 = z1 - dt * z1 / tau_z1 + noise_z1
    :z2 = z2 - dt * z2 / tau_z2 + noise_z2
    :z3 = z3 - dt * z3 / tau_z3 + noise_z3
    :z4 = z4 - dt * z4 / tau_z4 + noise_z4
    :z5 = z5 - dt * z5 / tau_z5 + noise_z5
    :z6 = z6 - dt * z6 / tau_z6 + noise_z6
    :z7 = z7 - dt * z7 / tau_z7 + noise_z7
    
    VERBATIM
    return 0;
    ENDVERBATIM
} : end PROCEDURE states()


PROCEDURE rates(Vm (mV)) { 
    LOCAL a,b,m3_inf,n4_inf,q10,sum,one_minus_m,one_minus_h,one_minus_n
    
    UNITSOFF
    
    q10 = 3.^((celsius-6.3)/10.)
    
    ::: SODIUM :::
    
    : 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
    
    : Exact
    mu_z1 = exp(-dt/tau_z1)
    mu_z2 = exp(-dt/tau_z2)
    mu_z3 = exp(-dt/tau_z3)
    mu_z4 = exp(-dt/tau_z4)
    mu_z5 = exp(-dt/tau_z5)
    mu_z6 = exp(-dt/tau_z6)
    mu_z7 = exp(-dt/tau_z7)
    noise_z1 = sqrt(var_z1 * (1-mu_z1*mu_z1)) * normrand(0,1)
    noise_z2 = sqrt(var_z2 * (1-mu_z2*mu_z2)) * normrand(0,1)
    noise_z3 = sqrt(var_z3 * (1-mu_z3*mu_z3)) * normrand(0,1)
    noise_z4 = sqrt(var_z4 * (1-mu_z4*mu_z4)) * normrand(0,1)
    noise_z5 = sqrt(var_z5 * (1-mu_z5*mu_z5)) * normrand(0,1)
    noise_z6 = sqrt(var_z6 * (1-mu_z6*mu_z6)) * normrand(0,1)
    noise_z7 = sqrt(var_z7 * (1-mu_z7*mu_z7)) * normrand(0,1)
    : Euler-Maruyama
    :noise_z1 = sqrt(2 * dt * var_z1 / tau_z1) * normrand(0,1)
    :noise_z2 = sqrt(2 * dt * var_z2 / tau_z2) * normrand(0,1)
    :noise_z3 = sqrt(2 * dt * var_z3 / tau_z3) * normrand(0,1)
    :noise_z4 = sqrt(2 * dt * var_z4 / tau_z4) * normrand(0,1)
    :noise_z5 = sqrt(2 * dt * var_z5 / tau_z5) * normrand(0,1)
    :noise_z6 = sqrt(2 * dt * var_z6 / tau_z6) * normrand(0,1)
    :noise_z7 = sqrt(2 * dt * var_z7 / tau_z7) * normrand(0,1)
    
    ::: POTASSIUM :::
    
    : 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
    tau_y1 = tau_n
    tau_y2 = tau_n/2
    tau_y3 = tau_n/3
    tau_y4 = tau_n/4
    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
    : Exact
    mu_y1 = exp(-dt/tau_y1)
    mu_y2 = exp(-dt/tau_y2)
    mu_y3 = exp(-dt/tau_y3)
    mu_y4 = exp(-dt/tau_y4)
    noise_y1 = sqrt(var_y1 * (1-mu_y1*mu_y1)) * normrand(0,1)
    noise_y2 = sqrt(var_y2 * (1-mu_y2*mu_y2)) * normrand(0,1)
    noise_y3 = sqrt(var_y3 * (1-mu_y3*mu_y3)) * normrand(0,1)
    noise_y3 = sqrt(var_y3 * (1-mu_y3*mu_y3)) * normrand(0,1)
    : Euler-Maruyama
    :noise_y1 = sqrt(2 * dt * var_y1 / tau_y1) * normrand(0,1)
    :noise_y2 = sqrt(2 * dt * var_y2 / tau_y2) * normrand(0,1)
    :noise_y3 = sqrt(2 * dt * var_y3 / tau_y3) * normrand(0,1)
    :noise_y4 = sqrt(2 * dt * var_y4 / tau_y4) * normrand(0,1)
    
    UNITSON
}

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
}


Loading data, please wait...