Fast sodium channel gating in mossy fiber axons (Schmidt-Hieber et al. 2010)

 Download zip file   Auto-launch 
Help downloading and running models
Accession:128079
"... To study the mechanisms underlying AP initiation in unmyelinated hippocampal mossy fibers of adult mice, we recorded sodium currents in axonal and somatic membrane patches. We demonstrate that sodium channel density in the proximal axon is ~5 times higher than in the soma. Furthermore, sodium channel activation and inactivation are ~2 times faster. Modeling revealed that the fast activation localized the initiation site to the proximal axon even upon strong synaptic stimulation, while fast inactivation contributed to energy-efficient membrane charging during APs. ..."
Reference:
1 . Schmidt-Hieber C, Bischofberger J (2010) Fast sodium channel gating supports localized and efficient axonal action potential initiation. J Neurosci 30:10233-42 [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:
Cell Type(s): Dentate gyrus granule cell;
Channel(s):
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: NEURON;
Model Concept(s): Action Potential Initiation;
Implementer(s):
Search NeuronDB for information about:  Dentate gyrus granule cell;
/
na8st
dat
hoc
mod
README.txt
ap.py
ap_utils.py
mech.py
qmatrix.py
stfio_plot.py
                            
import numpy as np

from scipy import fftpack
from scipy.interpolate import fitpack

import sys

import ap

class RunSettings(object):
    def __init__(self, gating=None, uniform_g=False, uniform_kin=False,
                 has_donnan=True,
                 gna_soma = 18.8,
                 gna_prox_axon = 94.0,
                 gna_distal_axon = 38.6152,
                 gk_scale_axon = 3.0, gk_scale_soma = 1.0, 
                 gk_axon = 0.004, gk_soma = 0.004, 
                 gk_distal_axon = 0.010):
        # Nav gating model
        self.gating = gating
        # Uniform gNa density
        self.uniform_g = uniform_g
        # Uniform gNa kinetics
        self.uniform_kin = uniform_kin
        # Apply correction for Donnan potential
        self.has_donnan = has_donnan
        # Somatic gNa
        self.gna_soma = gna_soma
        # Mean gNa in proximal axon (up to 40 um from soma)
        self.gna_prox_axon = gna_prox_axon
        # Distal axon gNa:
        self.gna_distal_axon = gna_distal_axon
        # Scaling factor for axonal Kv activation rates:
        self.gk_scale_axon = gk_scale_axon
        # Scaling factor for somatic Kv activation rates:
        self.gk_scale_soma = gk_scale_soma
        # Axonal Kv conductance density:
        self.gk_axon = gk_axon
        # Somatic Kv conductance density:
        self.gk_soma = gk_soma
        # Distal axonal Kv conductance density:
        self.gk_distal_axon = gk_distal_axon 

def fgaussColqu(x, f_c):
    """Eq. 5 from Colquhoun & Sigworth, p. 486 of the blue book"""
    # np.log(2.0)/2.0 = 0.34657359028
    return np.exp(-0.34657359028 * (x/f_c)**2)

def gaussian_filter(x, f_c, dt):
    """
    Gaussian filter according to Colquhoun & Sigworth (blue book).
    
    Arguments:
    x   -- Input data (1D numpy array)
    f_c -- Cutoff frequency in kHz (-3 dB)
    dt  -- sampling interval in ms

    Returns:
    x convolved with a Gaussian filter kernel.
    """

    xf = fftpack.rfft(x)
    # returns an array of frequencies corresponding to the indices of x:
    f = fftpack.rfftfreq(len(x), dt) 
    xf *= fgaussColqu(f, f_c)
    return fftpack.irfft(xf)

def rec_parent(ret_list, cell):
    sr = ap.h.SectionRef(sec = ret_list[-1])
    for candidate in cell.axo:
        if ap.h.secname(sec=sr.parent) == ap.h.secname(sec=candidate):
            ret_list.append(candidate)
            rec_parent(ret_list, cell)

def find_longest_axon(cell):
    max_dist = 0
    # set distance origin to border of soma:
    ap.h.distance(0,0.0, sec=cell.somaBorderLoc.secRef.sec)
    for axon in cell.axo:
        end_dist = ap.h.distance(1.0, sec=axon)
        if end_dist > max_dist:
            term_axon = axon
            max_dist = end_dist
    # recursively find parents:
    ret_list = [ term_axon, ]
    rec_parent(ret_list, cell)
    return ret_list

def maxRise(ap_wave, left, right):
    """
    returns the maximal slope of rise within the vector $o1 between indices $2 and $3
    (typically: beginning of event ($2) to index of peak ($3))
    adopted from stimfit
    """

    # Maximal rise
    maxRise=abs(ap_wave[right]-ap_wave[right-1])
    i=right-1
    while True:
        diff=abs(ap_wave[i]-ap_wave[i-1])
        if (maxRise < diff):
            maxRise=diff
        i -= 1
       	if (i<=left): break
    return maxRise

def maxDecay(ap_wave, left, right): 
    """
    returns the maximal slope of decay within the vector $o1 between indices $2 and $3
    (typically: index of peak ($2) to end of event ($3))
    adopted from stimfit
    """
    # {local left,right,maxDecay,i,diff
    if (left<0): left=0
    if (left > len(ap_wave)-3): left = len(ap_wave)-3
    if (right<0): right=0
    if (right > len(ap_wave)-1): right=len(ap_wave)-1
    # Maximal decay
    maxDecay=abs(ap_wave[left+1]-ap_wave[left])
    i=left+2
    while True:
        diff=abs(ap_wave[i]-ap_wave[i-1])
        if (maxDecay<diff):
            maxDecay=diff
        i+=1
        if (i>=right): break
    return maxDecay

def t50(ap_wave, peak_index, base, peak):
    """
    Returns the full width at half-maximal amplitude (FWHM) of an event in vector $o1.
    The index of the (estimated) peak of the event should be given as $2, baseline as $3, peak as $4
    adopted from stimfit
    """
    center = peak_index
    ampl = peak-base

    # walk left from center until HMA is reached:
    t50LeftId=center
    while True:
        t50LeftId-=1
        if t50LeftId <= 0:
            break
        if abs(ap_wave[t50LeftId]-base) < abs(0.5 * ampl):
            break

    # walk right:
    t50RightId=center
    while True:
        t50RightId += 1
        if t50RightId >= len(ap_wave):
            break
        if abs(ap_wave[t50RightId]-base)<abs(0.5 * ampl):
            break

    if t50LeftId >= len(ap_wave)-1:
        return 0
    if t50RightId < 0 or t50RightId >= len(ap_wave):
        return 0
    # calculation of real values by linear interpolation: 
    # Left side
    yLong2=ap_wave[t50LeftId+1]
    yLong1=ap_wave[t50LeftId]

    if yLong2-yLong1 != 0:
        t50LeftReal=(t50LeftId+abs((0.5*ampl-(yLong1-base))/(yLong2-yLong1)))
    else:
        t50LeftReal=t50LeftId

    # Right side
    yLong2=ap_wave[t50RightId]
    yLong1=ap_wave[t50RightId-1]

    if yLong2-yLong1 != 0:
        t50RightReal=t50RightId-abs((0.5*ampl-(yLong2-base))/abs(yLong2-yLong1))
    else:
        t50RightReal=t50RightId

    return t50RightReal-t50LeftReal

def analyse_ap(ap_wave, dt, f_c=0, resample=1):
    """
    Analyse an AP; return maximal rise, decay as a tuple
    f_c:      cutoff frequency for Gaussian filter 
              f_c == 0 means no filtering
    resample: factor for new sampling interval for resampling; i.e. dt_new = dt*resampling
              uses B-spline interpolation at present
              resampling==1 means no resampling
    """
    # local rise_axon, rise_soma, decay_axon, decay_soma, peak_index_axon, peak_index_soma, error, endInj, startInj
    if f_c > 0:
        pad1 = np.array([ap_wave[0] for i in range(1000)])
        pad2 = np.array([ap_wave[-1] for i in range(1000)])
        filtered = gaussian_filter(np.concatenate((pad1, ap_wave, pad2)), f_c, dt)[1000:1000+len(ap_wave)]
    else:
        filtered = ap_wave.copy()

    if resample != 1:
        # flin = interpolate.interp1d(np.arange(0,len(filtered),dtype=np.double)*dt, filtered)
        fspline = fitpack.splrep(np.arange(0,len(filtered),dtype=np.double)*dt, filtered)
        dt_ip = dt*resample
        y_ip = fitpack.splev(np.arange(0,len(filtered)/resample)*dt_ip, fspline)
    else:
        y_ip = filtered.copy()
        dt_ip = dt

    peak_index= y_ip.argmax()
    rise = maxRise(y_ip, 1, peak_index)/dt_ip
    decay = maxDecay(y_ip, peak_index, len(y_ip)-2)/dt_ip
    fwhm = t50(y_ip, peak_index, y_ip[0], y_ip.max()) * dt_ip
    return (rise, decay, fwhm, y_ip.max()-y_ip[0])

def whereis(wave, value):
    """
    returns the interpolated index of a vector ($o1) where
    $2 is found for the first time
    """
    # {local n,retIndex,fromtop,frombottom,m,c,x0,x1,y0,y1
    retIndex=0
    fromtop=False

    # coming from top or bottom?
    if (wave[0] > value):
	fromtop=True
    for n in range(len(wave)):
	if fromtop:
	    if (wave[n]<value):
		retIndex=n
		break
        elif (wave[n]>value):
            retIndex=n
            break
    if retIndex==0:
	sys.stderr.write("Value not found in whereis()\n")
	return 0

    # linear interpolation:
    x0=retIndex-1
    x1=retIndex
    y0=wave[x0]
    y1=wave[x1]
    m=(y1-y0)/(x1-x0)
    c=y0-m*x0
    return (value-c)/m

def init_h(bleb=False, dt=0.005, axon_seg=11, axon_lim=6, T=24.0, v_init=-80.0):
    """
    Initialize NEURON

    Arguments:
    bleb --     If True, cut the axon and add a spherical bleb
    dt --       Integration time step
    axon_seg -- Multiplication factor to increase nseg in the 
                proximal axon
    axon_lim -- Axonal section index up to which axon_seg is applied

    Returns:
    Cell object
    """

    ap.h("""load_file("stdrun.hoc")
            load_file("./hoc/config.hoc")
         """)
    
    ap.h.celsius = T
    ap.h.dt = dt
    ap.h.steps_per_ms = 1.0/ap.h.dt
    ap.h.v_init = v_init
    
    ap.h("""load_file("./hoc/membrane.hoc")""")

    arg = int(bleb)
    cell = ap.h.cell_10(arg) # 0 = no bleb, 1 = with bleb

    # Increase nseg in proximal axon:
    for (n_a,axon) in enumerate(cell.axo):
        if (n_a < axon_lim):
            axon.nseg *= axon_seg

    return cell

def fsigm(x, amp, center, slope):
    """
    Sigmoidal function
    """
    
    return amp-amp/(1.0+np.exp((x-center)/slope))

Loading data, please wait...