Neuron Hoc Tutorial

Example of how to use BMTOOL single cell analysis with neuron hoc files¶

By Gregory Glickert

First we make sure our modfiles are compiled

In [1]:

import os

# if already compiled then lets delete the folder and force a recompile
if os.path.isdir('modfiles/x86_64'):
    os.system("rm -rf modfiles/x86_64 ")
# compile the mod files
if not os.path.isdir("modfiles/x86_64"):
    os.chdir('modfiles')
    os.system("nrnivmodl")
    os.chdir("..")

import os # if already compiled then lets delete the folder and force a recompile if os.path.isdir('modfiles/x86_64'): os.system("rm -rf modfiles/x86_64 ") # compile the mod files if not os.path.isdir("modfiles/x86_64"): os.chdir('modfiles') os.system("nrnivmodl") os.chdir("..")

/home/gjgpb9/miniconda3/envs/bmtk/bin/nrnivmodl:10: DeprecationWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html
  from pkg_resources import working_set

/home/gjgpb9/cortex_modeling/bmtool/examples/single_cell/Neuron_hoc/modfiles
Mod files: "./AMPA_NMDA_STP.mod" "./cadad.mod" "./cal2.mod" "./can_mig.mod" "./exp2syn_stp.mod" "./GABA_A_STP.mod" "./gap.mod" "./Gfluct.mod" "./h_kole.mod" "./imCA3.mod" "./kap_BS.mod" "./kBK.mod" "./kdmc_BS.mod" "./kdr_BS.mod" "./kdrCA3.mod" "./kdrinter.mod" "./leak.mod" "./nainter.mod" "./napCA3.mod" "./natCA3.mod" "./nax_BS.mod" "./vecevent_coreneuron.mod"

Creating 'x86_64' directory for .o files.

 -> Compiling mod_func.cpp
 -> NMODL ../AMPA_NMDA_STP.mod
 -> NMODL ../cadad.mod
 -> NMODL ../cal2.mod
 -> NMODL ../can_mig.mod
Warning: Default 2 of PARAMETER cao will be ignored and set by NEURON.
Warning: Default 5e-05 of PARAMETER cai will be ignored and set by NEURON.
 -> NMODL ../exp2syn_stp.mod
 -> NMODL ../GABA_A_STP.mod
 -> NMODL ../gap.mod
Warning: Default 2 of PARAMETER cao will be ignored and set by NEURON.
Warning: Default 5e-05 of PARAMETER cai will be ignored and set by NEURON.
 -> NMODL ../Gfluct.mod
 -> NMODL ../h_kole.mod
 -> NMODL ../imCA3.mod
 -> NMODL ../kap_BS.mod
 -> NMODL ../kBK.mod
 -> NMODL ../kdmc_BS.mod
 -> NMODL ../kdr_BS.mod
 -> NMODL ../kdrCA3.mod
 -> NMODL ../kdrinter.mod
 -> NMODL ../leak.mod
 -> NMODL ../nainter.mod
 -> NMODL ../napCA3.mod
Warning: Default -80 of PARAMETER ek will be ignored and set by NEURON.
 -> NMODL ../natCA3.mod
Warning: Default 45 of PARAMETER ena will be ignored and set by NEURON.
 -> NMODL ../nax_BS.mod
 -> NMODL ../vecevent_coreneuron.mod
 -> Compiling AMPA_NMDA_STP.c
 -> Compiling cadad.c
Notice: ARTIFICIAL_CELL is a synonym for POINT_PROCESS which hints that it
only affects and is affected by discrete events. As such it is not
located in a section and is not associated with an integrator
 -> Compiling cal2.c
 -> Compiling can_mig.c
 -> Compiling exp2syn_stp.c
 -> Compiling GABA_A_STP.c
 -> Compiling gap.c
 -> Compiling Gfluct.c
 -> Compiling h_kole.c
 -> Compiling imCA3.c
 -> Compiling kap_BS.c
 -> Compiling kBK.c
 -> Compiling kdmc_BS.c
 -> Compiling kdr_BS.c
 -> Compiling kdrCA3.c
 -> Compiling kdrinter.c
 -> Compiling leak.c
 -> Compiling nainter.c
 -> Compiling napCA3.c
 -> Compiling natCA3.c
 -> Compiling nax_BS.c
 -> Compiling vecevent_coreneuron.c

Translating AMPA_NMDA_STP.mod into /home/gjgpb9/cortex_modeling/bmtool/examples/single_cell/Neuron_hoc/modfiles/x86_64/AMPA_NMDA_STP.c
Translating cal2.mod into /home/gjgpb9/cortex_modeling/bmtool/examples/single_cell/Neuron_hoc/modfiles/x86_64/cal2.c
Translating cadad.mod into /home/gjgpb9/cortex_modeling/bmtool/examples/single_cell/Neuron_hoc/modfiles/x86_64/cadad.c
Thread Safe
Thread Safe
Thread Safe
Translating can_mig.mod into /home/gjgpb9/cortex_modeling/bmtool/examples/single_cell/Neuron_hoc/modfiles/x86_64/can_mig.c
Thread Safe
Translating exp2syn_stp.mod into /home/gjgpb9/cortex_modeling/bmtool/examples/single_cell/Neuron_hoc/modfiles/x86_64/exp2syn_stp.c
Translating GABA_A_STP.mod into /home/gjgpb9/cortex_modeling/bmtool/examples/single_cell/Neuron_hoc/modfiles/x86_64/GABA_A_STP.c
Translating gap.mod into /home/gjgpb9/cortex_modeling/bmtool/examples/single_cell/Neuron_hoc/modfiles/x86_64/gap.c
Thread Safe
Thread Safe
Thread Safe
Translating Gfluct.mod into /home/gjgpb9/cortex_modeling/bmtool/examples/single_cell/Neuron_hoc/modfiles/x86_64/Gfluct.c
Notice: This mechanism cannot be used with CVODE
Thread Safe
Translating h_kole.mod into /home/gjgpb9/cortex_modeling/bmtool/examples/single_cell/Neuron_hoc/modfiles/x86_64/h_kole.c
Translating imCA3.mod into /home/gjgpb9/cortex_modeling/bmtool/examples/single_cell/Neuron_hoc/modfiles/x86_64/imCA3.c
Translating kap_BS.mod into /home/gjgpb9/cortex_modeling/bmtool/examples/single_cell/Neuron_hoc/modfiles/x86_64/kap_BS.c
Thread Safe
Thread Safe
Thread Safe
Translating kBK.mod into /home/gjgpb9/cortex_modeling/bmtool/examples/single_cell/Neuron_hoc/modfiles/x86_64/kBK.c
Thread Safe
Translating kdr_BS.mod into /home/gjgpb9/cortex_modeling/bmtool/examples/single_cell/Neuron_hoc/modfiles/x86_64/kdr_BS.c
Translating kdmc_BS.mod into /home/gjgpb9/cortex_modeling/bmtool/examples/single_cell/Neuron_hoc/modfiles/x86_64/kdmc_BS.c
Translating kdrCA3.mod into /home/gjgpb9/cortex_modeling/bmtool/examples/single_cell/Neuron_hoc/modfiles/x86_64/kdrCA3.c
Thread Safe
Thread Safe
Thread Safe
Translating kdrinter.mod into /home/gjgpb9/cortex_modeling/bmtool/examples/single_cell/Neuron_hoc/modfiles/x86_64/kdrinter.c
Thread Safe
Translating leak.mod into /home/gjgpb9/cortex_modeling/bmtool/examples/single_cell/Neuron_hoc/modfiles/x86_64/leak.c
Translating nainter.mod into /home/gjgpb9/cortex_modeling/bmtool/examples/single_cell/Neuron_hoc/modfiles/x86_64/nainter.c
Translating napCA3.mod into /home/gjgpb9/cortex_modeling/bmtool/examples/single_cell/Neuron_hoc/modfiles/x86_64/napCA3.c
Thread Safe
Thread Safe
Thread Safe
Translating natCA3.mod into /home/gjgpb9/cortex_modeling/bmtool/examples/single_cell/Neuron_hoc/modfiles/x86_64/natCA3.c
Thread Safe
Translating nax_BS.mod into /home/gjgpb9/cortex_modeling/bmtool/examples/single_cell/Neuron_hoc/modfiles/x86_64/nax_BS.c
Translating vecevent_coreneuron.mod into /home/gjgpb9/cortex_modeling/bmtool/examples/single_cell/Neuron_hoc/modfiles/x86_64/vecevent_coreneuron.c
Thread Safe
Thread Safe

 => LINKING shared library ./libnrnmech.so
 => LINKING executable ./special LDFLAGS are:    -pthread
Successfully created x86_64/special

Init our profiler¶

we need to init our profiler this will make sure the templates and mechanisms are loaded

In [2]:

import numpy as np
from neuron import h
import matplotlib.pyplot as plt
from bmtool.singlecell import Profiler, Passive, CurrentClamp, FI, ZAP, run_and_plot
%matplotlib inline

template_dir = '.' #templates are in templates.hoc in this working dir
mechanism_dir = 'modfiles'

profiler = Profiler(template_dir=template_dir, mechanism_dir=mechanism_dir, dt = 0.05)

import numpy as np from neuron import h import matplotlib.pyplot as plt from bmtool.singlecell import Profiler, Passive, CurrentClamp, FI, ZAP, run_and_plot %matplotlib inline template_dir = '.' #templates are in templates.hoc in this working dir mechanism_dir = 'modfiles' profiler = Profiler(template_dir=template_dir, mechanism_dir=mechanism_dir, dt = 0.05)

Warning: no DISPLAY environment variable.
--No graphics will be displayed.

Then we will setup our settings for each function in the singlecell module. These setting may change depending on cell type and the biology data you have to constain off of. For now some basic ones are used for looking at principle cells¶

In [3]:

noise = False
post_init_function = 'insert_mechs(0)' if noise else None

basic_settings = {
    'Passive': {
        'celsius': 26.,
        'kwargs': {
            'inj_amp': -50.,
            'inj_delay': 1500.,
            'inj_dur': 1000.,
            'tstop': 2500.,
            'method': 'exp2' # can be exp2, exp or simple
        }
    },
    'CurrentClamp': {
        'celsius': 34.,
        'kwargs': {
            'post_init_function': post_init_function,
            'inj_amp': 350.,
            'inj_delay': 1500.,
            'inj_dur': 1000.,
            'tstop': 3000.,
            'threshold': 0.
        }
    },
    'ZAP': {
        'celsius': 34.,
        'kwargs': {
            'post_init_function': post_init_function,
            'inj_amp': 100.,
            'inj_delay': 1000.,
            'inj_dur': 15000.,
            'tstop': 15500.,
            'fstart': 0.,
            'fend': 15.,
            'chirp_type': 'linear'
        }
    },
    'FI': {
        'celsius': 34.,
        'kwargs': {
            'post_init_function': post_init_function,
            'i_start': 0.,
            'i_stop': 2000.,
            'i_increment': 100.,
            'tstart': 1500.
        }
    }
}

Cell = 'CP_Cell'

noise = False post_init_function = 'insert_mechs(0)' if noise else None basic_settings = { 'Passive': { 'celsius': 26., 'kwargs': { 'inj_amp': -50., 'inj_delay': 1500., 'inj_dur': 1000., 'tstop': 2500., 'method': 'exp2' # can be exp2, exp or simple } }, 'CurrentClamp': { 'celsius': 34., 'kwargs': { 'post_init_function': post_init_function, 'inj_amp': 350., 'inj_delay': 1500., 'inj_dur': 1000., 'tstop': 3000., 'threshold': 0. } }, 'ZAP': { 'celsius': 34., 'kwargs': { 'post_init_function': post_init_function, 'inj_amp': 100., 'inj_delay': 1000., 'inj_dur': 15000., 'tstop': 15500., 'fstart': 0., 'fend': 15., 'chirp_type': 'linear' } }, 'FI': { 'celsius': 34., 'kwargs': { 'post_init_function': post_init_function, 'i_start': 0., 'i_stop': 2000., 'i_increment': 100., 'tstart': 1500. } } } Cell = 'CP_Cell'

Passive Properties¶

In [4]:

proc = basic_settings['Passive']
h.celsius = proc['celsius']
sim = Passive(Cell, **proc['kwargs'])
title = 'Passive Cell Current Injection'
xlabel = 'Time (ms)'
ylabel = 'Membrane Potential (mV)'

proc = basic_settings['Passive'] h.celsius = proc['celsius'] sim = Passive(Cell, **proc['kwargs']) title = 'Passive Cell Current Injection' xlabel = 'Time (ms)' ylabel = 'Membrane Potential (mV)'

Injection location: CP_Cell[0].soma[0](0.5)
Recording: CP_Cell[0].soma[0](0.5)._ref_v

we will have different plots for the different methods

In [5]:

if sim.method =='exp2':
    X, Y = run_and_plot(sim, title, xlabel, ylabel, plot_injection_only=True)
    plt.gca().plot(*sim.double_exponential_fit(), 'r:', label='double exponential fit')
    plt.legend()
elif sim.method =='exp':
    X, Y = run_and_plot(sim, title, xlabel, ylabel, plot_injection_only=True)
    plt.gca().plot(*sim.single_exponential_fit(), 'r:', label='single exponential fit')
    plt.legend()
else:
    X, Y = run_and_plot(sim, title, xlabel, ylabel, plot_injection_only=True)
plt.show()

if sim.method =='exp2': X, Y = run_and_plot(sim, title, xlabel, ylabel, plot_injection_only=True) plt.gca().plot(*sim.double_exponential_fit(), 'r:', label='double exponential fit') plt.legend() elif sim.method =='exp': X, Y = run_and_plot(sim, title, xlabel, ylabel, plot_injection_only=True) plt.gca().plot(*sim.single_exponential_fit(), 'r:', label='single exponential fit') plt.legend() else: X, Y = run_and_plot(sim, title, xlabel, ylabel, plot_injection_only=True) plt.show()

Running simulation for passive properties...

V Rest: -67.74 (mV)
Resistance: 86.03 (MOhms)
Membrane time constant: 39.95 (ms)

V_rest Calculation: Voltage taken at time 1500.0 (ms) is
-67.74 (mV)

R_in Calculation: dV/dI = (v_final-v_rest)/(i_final-i_start)
(-72.04 - (-67.74)) / (-0.05 - 0)
4.30 (mV) / 0.05 (nA) = 86.03 (MOhms)

Tau Calculation: Fit a double exponential curve to the membrane potential response
f(t) = a0 + a1*exp(-t/tau1) + a2*exp(-t/tau2)
Constrained by initial value: f(0) = a0 + a1 + a2 = v_rest
Fit parameters: (a0, a1, a2, tau1, tau2) = (-72.03, -0.48, 4.77, 39.95, 11.90)
Membrane time constant is determined from the slowest exponential term: 39.95 (ms)

Sag potential: v_sag = v_peak - v_final = -0.08 (mV)
Normalized sag potential: v_sag / (v_peak - v_rest) = 0.019

Current Injection¶

In [6]:

proc = basic_settings['CurrentClamp']
h.celsius = proc['celsius']
sim2 = CurrentClamp(Cell, **proc['kwargs'])
title = 'Current Injection'
xlabel = 'Time (ms)'
ylabel = 'Membrane Potential (mV)'

proc = basic_settings['CurrentClamp'] h.celsius = proc['celsius'] sim2 = CurrentClamp(Cell, **proc['kwargs']) title = 'Current Injection' xlabel = 'Time (ms)' ylabel = 'Membrane Potential (mV)'

Injection location: CP_Cell[1].soma[0](0.5)
Recording: CP_Cell[1].soma[0](0.5)._ref_v

In [7]:

X, Y = run_and_plot(sim2, title, xlabel, ylabel, plot_injection_only=True)
plt.show()

X, Y = run_and_plot(sim2, title, xlabel, ylabel, plot_injection_only=True) plt.show()

Current clamp simulation running...

Number of spikes: 27

Impedance Amplitude Profile (ZAP)¶

In [8]:

proc = basic_settings['ZAP']
h.celsius = proc['celsius']
sim3 = ZAP(Cell, **proc['kwargs'])
title = 'ZAP Response'
xlabel = 'Time (ms)'
ylabel = 'Membrane Potential (mV)'

proc = basic_settings['ZAP'] h.celsius = proc['celsius'] sim3 = ZAP(Cell, **proc['kwargs']) title = 'ZAP Response' xlabel = 'Time (ms)' ylabel = 'Membrane Potential (mV)'

Injection location: CP_Cell[2].soma[0](0.5)
Recording: CP_Cell[2].soma[0](0.5)._ref_v

In [9]:

X, Y = run_and_plot(sim3, title, xlabel, ylabel, plot_injection_only=True)

plt.figure()
plt.plot(X, sim3.zap_vec)
plt.title('ZAP Current')
plt.xlabel('Time (ms)')
plt.ylabel('Current Injection (nA)')

plt.figure()
plt.plot(*sim3.get_impedance(smooth=9))
plt.title('Impedance Amplitude Profile')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Impedance (MOhms)')
plt.show()

X, Y = run_and_plot(sim3, title, xlabel, ylabel, plot_injection_only=True) plt.figure() plt.plot(X, sim3.zap_vec) plt.title('ZAP Current') plt.xlabel('Time (ms)') plt.ylabel('Current Injection (nA)') plt.figure() plt.plot(*sim3.get_impedance(smooth=9)) plt.title('Impedance Amplitude Profile') plt.xlabel('Frequency (Hz)') plt.ylabel('Impedance (MOhms)') plt.show()

ZAP current simulation running...

Chirp current injection with frequency changing from 0 to 15 Hz over 15 seconds
Impedance is calculated as the ratio of FFT amplitude of membrane voltage to FFT amplitude of chirp current

Resonant Peak Frequency: 0.333 (Hz)

F-I Curve¶

In [10]:

proc = basic_settings['FI']
h.celsius = proc['celsius']
sim4 = FI(Cell, **proc['kwargs'])
title = 'FI Curve'
xlabel = 'Injection (nA)'
ylabel = '# Spikes'

proc = basic_settings['FI'] h.celsius = proc['celsius'] sim4 = FI(Cell, **proc['kwargs']) title = 'FI Curve' xlabel = 'Injection (nA)' ylabel = '# Spikes'

Injection location: CP_Cell[22].soma[0](0.5)
Recording: CP_Cell[22].soma[0](0.5)._ref_v

In [11]:

X, Y = run_and_plot(sim4, title, xlabel, ylabel)
plt.show()

X, Y = run_and_plot(sim4, title, xlabel, ylabel) plt.show()

Running simulations for FI curve...

Results
    Injection (nA):  number of spikes
0               0.0                 0
1               0.1                 0
2               0.2                15
3               0.3                23
4               0.4                31
5               0.5                40
6               0.6                49
7               0.7                59
8               0.8                70
9               0.9                 4
10              1.0                 2
11              1.1                 2
12              1.2                 2
13              1.3                 2
14              1.4                 2
15              1.5                 1
16              1.6                 1
17              1.7                 1
18              1.8                 1
19              1.9                 1

You can then take the output of the FI curve and find some simple properties of it

In [12]:

def find_slope(X, Y):
    """Finds the slope and intercept using least squares
    
    """
    # Create the A matrix with a column of X and a column of ones (for the intercept term)
    A = np.vstack([X, np.ones_like(X)]).T

    # Perform least-squares fitting to find the slope (m) and intercept (c)
    m, c = np.linalg.lstsq(A, Y, rcond=None)[0]

    # Return the slope (m) and intercept (c)
    return m, c

def fit_slope_to_max_y(X, Y):
    """Fits a line between Y > 0 and the point where Y reaches its maximum"""
    # Find where Y > 0
    positive_data = Y > 0
    X_positive = X[positive_data]
    Y_positive = Y[positive_data]

    # Find the index of the maximum Y value
    max_idx = np.argmax(Y_positive)

    # Only use data from where Y is increasing up to the max value
    X_fit = X_positive[:max_idx + 1]
    Y_fit = Y_positive[:max_idx + 1]

    # Find the slope and intercept for this range
    m, c = find_slope(X_fit, Y_fit)

    return m, c, X_fit, Y_fit

def find_rheobase(X,Y):
    """Rheobase is min current to cause a spike"""
    non_zero_indices = np.nonzero(Y)[0]
    rheobase = X[non_zero_indices[0]]
    print("Rheobase: " + str(rheobase) + 'pA')

X = X * 1000 #nA to pA
m, c, X_fit, Y_fit = fit_slope_to_max_y(X, Y)
plt.scatter(X,Y)
plt.plot(X_fit, m*X_fit + c, label=f'Fit Line: y = {m:.2f}x + {c:.2f}', color='red')
plt.title("Fitted slope to FI curve")
plt.legend()
plt.show()

find_rheobase(X,Y)

def find_slope(X, Y): """Finds the slope and intercept using least squares """ # Create the A matrix with a column of X and a column of ones (for the intercept term) A = np.vstack([X, np.ones_like(X)]).T # Perform least-squares fitting to find the slope (m) and intercept (c) m, c = np.linalg.lstsq(A, Y, rcond=None)[0] # Return the slope (m) and intercept (c) return m, c def fit_slope_to_max_y(X, Y): """Fits a line between Y > 0 and the point where Y reaches its maximum""" # Find where Y > 0 positive_data = Y > 0 X_positive = X[positive_data] Y_positive = Y[positive_data] # Find the index of the maximum Y value max_idx = np.argmax(Y_positive) # Only use data from where Y is increasing up to the max value X_fit = X_positive[:max_idx + 1] Y_fit = Y_positive[:max_idx + 1] # Find the slope and intercept for this range m, c = find_slope(X_fit, Y_fit) return m, c, X_fit, Y_fit def find_rheobase(X,Y): """Rheobase is min current to cause a spike""" non_zero_indices = np.nonzero(Y)[0] rheobase = X[non_zero_indices[0]] print("Rheobase: " + str(rheobase) + 'pA') X = X * 1000 #nA to pA m, c, X_fit, Y_fit = fit_slope_to_max_y(X, Y) plt.scatter(X,Y) plt.plot(X_fit, m*X_fit + c, label=f'Fit Line: y = {m:.2f}x + {c:.2f}', color='red') plt.title("Fitted slope to FI curve") plt.legend() plt.show() find_rheobase(X,Y)

Slope: 0.5263157894736842
Rheobase: 0.2nA

Download NotebookDownload Folder