Channel models (Notebook)

In this tutorial, we will learn how to work with ion channel models.

Setup

!pip install dendrotweaks --quiet

If you are using Google Colab, you might also need to restart the session as the installation downgraded some packages (numpy). You can do it manually or programmatically as shown below:

# import os
# os.kill(os.getpid(), 9)

import numpy as np
import matplotlib.pyplot as plt
import os

Let’s begin by importing the standard libraries and the dendrotweaks library.

import dendrotweaks as dd
dd.__version__

'0.5.1'

dd.apply_dark_theme() # Set the theme for the plots

os.makedirs('examples', exist_ok=True)
if not os.listdir('examples/'):
    print("Downloading example data...")
    dd.download_example_data('examples')

Converting MOD to Python

The MODFileConverter class conveniently encapsulates the process of converting MOD files into Python code. Let’s create an instance of this class.

converter = dd.biophys.io.MODFileConverter()

converter.reader, converter.parser, converter.generator

(<dendrotweaks.biophys.io.reader.MODFileReader at 0x7f458a79e210>,
 <dendrotweaks.biophys.io.parser.MODFileParser at 0x7f455e6c9090>,
 <dendrotweaks.biophys.io.code_generators.PythonCodeGenerator at 0x7f455e6c9010>)

path_to_mod_file = os.path.join('examples', 'Toy', 'biophys', 'mod', 'Kv.mod')
path_to_python = os.path.join('examples', 'Toy','biophys', 'python', 'Kv.py')
path_to_template = os.path.join('examples', 'Templates', 'default.py')

converter.convert(path_to_mod_file, path_to_python, path_to_template)

Warning: Procedure compute_kinetic_variables has no parameters! Expected 'v' or 'cai'. Defaulting to 'v'.
Saved content to examples/Toy/biophys/python/Kv.py

The conversion process involves reading a MOD file, parsing its content to build an abstract syntax tree representation, and then generating new code using a Jinja template.

The MODReader class reads the files, performs basic preprocessing steps like removing comments and empty lines, and splits the file content into distinct blocks.

blocks = converter.reader.get_blocks()

Split content into blocks:
    1 - TITLE
    1 - COMMENT
    1 - NEURON
    1 - UNITS
    1 - PARAMETER
    1 - ASSIGNED
    1 - STATE
    1 - BREAKPOINT
    1 - DERIVATIVE
    1 - INITIAL
    1 - FUNCTION
    1 - PROCEDURE
    0 - KINETIC

blocks

{'TITLE': ['TITLE Kv_Park_ref\n'],
 'COMMENT': ['COMMENT\n26 Ago 2002 Modification of original channel to allow variable time step and to correct an initialization error.\n    Done by Michael Hines(michael.hines@yale.e) and Ruggero Scorcioni(rscorcio@gmu.edu) at EU Advance Course in Computational Neuroscience. Obidos, Portugal\n\nkv.mod\n\nPotassium channel, Hodgkin-Huxley style kinetics\nKinetic rates based roughly on Sah et al. and Hamill et al. (1991)\n\nAuthor\nENDCOMMENT'],
 'NEURON': ['NEURON {\n\tSUFFIX Kv\n\tUSEION k READ ek WRITE ik\n\tRANGE gbar, i, v12, q\n}'],
 'UNITS': ['UNITS {\n\t(mA) = (milliamp)\n\t(mV) = (millivolt)\n\t(S)  = (siemens)\n\t(um) = (micron)\n}'],
 'PARAMETER': ['PARAMETER {\n\tgbar = 0.0   (S/cm2)                                                                                      \n\tRa   = 0.02  (/mV/ms) \n\tRb   = 0.006 (/mV/ms) \n\tv12  = 25    (mV)\t  \n\tq    = 9     (mV)     \n\ttemp = 23    (degC)   \n\tq10  = 2.3   (1)\t  \n}'],
 'ASSIGNED': ['ASSIGNED {\n\tv       (mV)\n\ti \t    (mA/cm2)\n\tik      (mA/cm2)\n\tgk      (S/cm2)\n\tek      (mV)\n\tninf    (1)\n\tntau    (ms)\n\ttadj    (1)\n\tcelsius (degC)\n}'],
 'STATE': ['STATE { n }'],
 'BREAKPOINT': ['BREAKPOINT {\n    SOLVE states METHOD cnexp\n\tgk = tadj * gbar * n\n\ti = gk * (v - ek)\n\tik = i\n}'],
 'DERIVATIVE': ["DERIVATIVE  states {   \n        rates()\n        n' =  (ninf - n)/ntau\n}"],
 'INITIAL': ['INITIAL { \n\trates()\n\tn = ninf\n}'],
 'FUNCTION': ['FUNCTION rateconst(v (mV), r (/mV/ms), th (mV), q (mV)) (/ms) {\n\trateconst = r * (v - th) / (1 - exp(-(v - th)/q))\n}'],
 'PROCEDURE': ['PROCEDURE rates() {\n\n\tLOCAL alpn, betn\n\n\ttadj = q10^((celsius - temp)/10(degC))\n\n\talpn = rateconst(v, Ra, v12, q)\n\tbetn = rateconst(v, -Rb, v12, -q)\n\n    ntau = 1 / (tadj * (alpn + betn))\n\tninf = alpn/(alpn + betn)\n}'],
 'KINETIC': []}

The MOD files are parsed using the PyParsing library. The grammar used for parsing can be found in the dendrotweaks/biophys/grammar.py file. The parser generates an Abstract Syntax Tree (AST) from the MOD file. The AST is a hierarchical representation of the MOD file that can be used to generate Python code.

ast = converter.parser.get_ast()

ast

{'TITLE': ' Kv_Park_ref',
 'COMMENT': '\n'
            '26 Ago 2002 Modification of original channel to allow variable '
            'time step and to correct an initialization error.\n'
            '    Done by Michael Hines(michael.hines@yale.e) and Ruggero '
            'Scorcioni(rscorcio@gmu.edu) at EU Advance Course in Computational '
            'Neuroscience. Obidos, Portugal\n'
            '\n'
            'kv.mod\n'
            '\n'
            'Potassium channel, Hodgkin-Huxley style kinetics\n'
            'Kinetic rates based roughly on Sah et al. and Hamill et al. '
            '(1991)\n'
            '\n'
            'Author\n',
 'NEURON': {'suffix': 'Kv',
            'useion': [{'ion': 'k', 'read': ['ek'], 'write': ['ik']}],
            'range': ['gbar', 'i', 'v12', 'q']},
 'UNITS': {'mA': 'milliamp', 'mV': 'millivolt', 'S': 'siemens', 'um': 'micron'},
 'PARAMETER': [{'name': 'gbar', 'value': 0.0, 'unit': 'S/cm2'},
               {'name': 'Ra', 'value': 0.02, 'unit': '/mV/ms'},
               {'name': 'Rb', 'value': 0.006, 'unit': '/mV/ms'},
               {'name': 'v12', 'value': 25, 'unit': 'mV'},
               {'name': 'q', 'value': 9, 'unit': 'mV'},
               {'name': 'temp', 'value': 23, 'unit': 'degC'},
               {'name': 'q10', 'value': 2.3, 'unit': '1'}],
 'ASSIGNED': [{'name': 'v', 'unit': 'mV'},
              {'name': 'i', 'unit': 'mA/cm2'},
              {'name': 'ik', 'unit': 'mA/cm2'},
              {'name': 'gk', 'unit': 'S/cm2'},
              {'name': 'ek', 'unit': 'mV'},
              {'name': 'ninf', 'unit': '1'},
              {'name': 'ntau', 'unit': 'ms'},
              {'name': 'tadj', 'unit': '1'},
              {'name': 'celsius', 'unit': 'degC'}],
 'STATE': {'n': {'power': 1}},
 'BREAKPOINT': {'solve_stmt': {'solve': 'states', 'method': 'cnexp'},
                'statements': [{'assigned_var': 'gk',
                                'expression': '(tadj * gbar) * n'},
                               {'assigned_var': 'i',
                                'expression': 'gk * (v - ek)'},
                               {'assigned_var': 'ik', 'expression': 'i'}]},
 'DERIVATIVE': {'name': 'states',
                'func_calls': [{'rates': []}],
                'statements': [{'assigned_var': 'n',
                                'expression': '(ninf - n) / ntau'}]},
 'INITIAL': {'statements': [{'rates': []},
                            {'assigned_var': 'n', 'expression': 'ninf'}]},
 'FUNCTION': [{'signature': {'name': 'rateconst',
                             'params': [{'name': 'v', 'unit': 'mV'},
                                        {'name': 'r', 'unit': '/mV/ms'},
                                        {'name': 'th', 'unit': 'mV'},
                                        {'name': 'q', 'unit': 'mV'}],
                             'returned_unit': '/ms'},
               'locals': [],
               'statements': [{'assigned_var': 'rateconst',
                               'expression': '(r * (v - th)) / (1 - '
                                             'np.exp((-(v - th) / q)))'}]}],
 'PROCEDURE': [{'signature': {'name': 'compute_kinetic_variables'},
                'locals': ['alpn', 'betn'],
                'statements': [{'assigned_var': 'tadj',
                                'expression': 'q10 ** ((celsius - temp) / 10)'},
                               {'assigned_var': 'alpn',
                                'expression': 'rateconst(v, Ra, v12, q)'},
                               {'assigned_var': 'betn',
                                'expression': 'rateconst(v, -Rb, v12, -q)'},
                               {'assigned_var': 'nTau',
                                'expression': '1 / (tadj * (alpn + betn))'},
                               {'assigned_var': 'nInf',
                                'expression': 'alpn / (alpn + betn)'}]}]}

The PythonCodeGenerator class generates Python code from the AST using a JINJA template file.

print(converter.generator.content)

# This Python channel class was automatically generated from a MOD file
# using DendroTweaks toolbox, dendrotweaks.dendrites.gr

import sys

from dendrotweaks.biophys.mechanisms import IonChannel
import numpy as np

class Kv(IonChannel):
    """
    Kv_Park_ref
    """

    def __init__(self, name="Kv"):
        super().__init__(name=name)
        self.params = {
            "gbar": 0.0,
            "Ra": 0.02,
            "Rb": 0.006,
            "v12": 25,
            "q": 9,
            "temp": 23,
            "q10": 2.3
            }
        self.range_params = {
            "gbar": 0.0,
            "v12": 25,
            "q": 9
            }
        self.states = {
            "n": 0.0
            }
        self._state_powers = {
            "n": {'power': 1}
            }
        self.ion = "k"
        self.current_name = "i_k"
        self.current_available = True
        self.independent_var_name = "v"
        self.temperature = 37

    def __getitem__(self, item):
        return self.params[item]

    def __setitem__(self, item, value):
        self.params[item] = value

    
    def compute_kinetic_variables(self, v):
        Ra = self.params["Ra"]
        Rb = self.params["Rb"]
        v12 = self.params["v12"]
        q = self.params["q"]
        
        alpn = self.rateconst(v, Ra, v12, q)
        betn = self.rateconst(v, -Rb, v12, -q)
        nTau = 1 / (self.tadj * (alpn + betn))
        nInf = alpn / (alpn + betn)
        return nInf, nTau
    
    
    def rateconst(self, v, r, th, q):
        
        rateconst = (r * (v - th)) / (1 - np.exp((-(v - th) / q)))
        return rateconst

Working with channels

We can create an independent channel object outside a model as follows:

from dendrotweaks.biophys.io import create_channel

path_to_mod_file = os.path.join('examples', 'Toy', 'biophys', 'mod', 'Kv.mod')
path_to_python = os.path.join('examples', 'Toy','biophys', 'python', 'Kv.py')
path_to_template = os.path.join('examples', 'Templates', 'default.py')

kv_channel = create_channel(path_to_mod_file, path_to_python, path_to_template)

Warning: Procedure compute_kinetic_variables has no parameters! Expected 'v' or 'cai'. Defaulting to 'v'.
Saved content to examples/Toy/biophys/python/Kv.py

fig, ax = plt.subplots(1, 2, figsize=(10, 5))
kv_channel.plot_kinetics(ax=ax)

Got data for v in range -100.0 to 100.0 at 37°C

However, typically we would first create a model and then add channels to it. Let’s create a simple neuron model.

path_to_model = os.path.join('examples', 'Toy')
print(f'Path to model: {path_to_model}')

Path to model: examples/Toy

model = dd.Model(path_to_model)

model.list_morphologies()

['simple-cell-ais', 'ball-and-stick', 'point', 'simple-cell']

model.load_morphology(file_name='simple-cell')

Sorted PointTree(root=Point(idx=0), num_nodes=17).
Apical dendrite is already aligned.
Extended 6 nodes.
Sorted PointTree(root=Point(idx=0), num_nodes=23).

model.sec_tree.plot(
    show_points=True,
    annotate=True,
)

On the model level, we can create an ion channel using the add_mechanism method. First, let’s add the default mechanisms to the model, which include the leak channel, calcium dynamics and synaptic mechanisms. Next, we add the sodium and potassium channels to the model.

model.add_default_mechanisms()
model.add_mechanism('Na')
model.add_mechanism('Kv')

Saved content to examples/Toy/biophys/python/Na.py
Warning: Procedure compute_kinetic_variables has no parameters! Expected 'v' or 'cai'. Defaulting to 'v'.
Saved content to examples/Toy/biophys/python/Kv.py

With these commands, we create Python objects from MOD files and add them to mechanisms. We also compile and load the MOD files, making them available in NEURON.

Once the ion channel is created, its kinetics can be visualized using the plot_kinetics method.

fig, ax = plt.subplots(1, 2, figsize=(10, 5))
model.mechanisms['Kv'].plot_kinetics(ax=ax)

Got data for v in range -100.0 to 100.0 at 37°C

Each mechanism object is an instance of the Mechanism class, which contains information about the mechanism, such as its name and parameters. To examine the parameters of a specific mechanism, we can use the params attribute:

model.mechanisms['Kv'].params

{'gbar': 0.0,
 'Ra': 0.02,
 'Rb': 0.006,
 'v12': 25,
 'q': 9,
 'temp': 23,
 'q10': 2.3}

After uploading the mechanisms in the previous step, we now need to insert them into specific domains. In this example, we insert each of the three available mechanisms to all domains.

for domain_name in model.domains:
    for mech_name in model.mechanisms:
        model.insert_mechanism(mech_name, domain_name, distribute=False)

At this point, we have inserted the mechanisms into the membrane and set the default parameters for the model. We can inspect the parameters of the inserted mechanisms using the df_params or params attribute.

model.df_params[model.df_params['Mechanism'] == 'Kv']

	Mechanism	Parameter	Group	Distribution	Distribution params
24	Kv	gbar_Kv	all	constant	{'value': 0.0}
25	Kv	v12_Kv	all	constant	{'value': 25}
26	Kv	q_Kv	all	constant	{'value': 9}

Let’s update channel conductances and run a simulation in out point neuron.

model.set_param('gbar_Leak', value=0.00015)
model.set_param('gbar_Na', value=0.03)
model.set_param('gbar_Na', group_name='somatic', value=0.05)
model.set_param('gbar_Kv', value=0.003)
model.set_param('gbar_Kv', group_name='somatic', value=0.005)

model.remove_all_stimuli()
model.remove_all_recordings()

soma = model.sec_tree[0]
model.add_recording(soma, loc=0.5)
model.add_iclamp(soma, loc=0.5, amp=0.1, delay=50, dur=200)

model.run(300)
fig, ax = plt.subplots(1, 1, figsize=(10, 4))
model.simulator.plot_voltage(ax=ax)
t, v_orig = model.simulator.t, model.simulator.recordings['v'][soma(0.5)]

Standardizing ion channel models

DendroTweaks provides a simple way to standardize voltage-gated ion channel models.

Assuming that we have an instance of a model with a custom sodium channel Na created before, we can standardize it with the model.standardize_channel method. The standard channel with the name stdNa (standard-Na) and very similar kinetics to the original channel model is created to replace the original channel model. This method also exports the standard channel model to a MOD file which is immediately loaded into the NEURON simulator.

fig, ax = plt.subplots(1, 2, figsize=(10, 5))
model.mechanisms['Na'].plot_kinetics(ax=ax)
model.standardize_channel('Na');
model.mechanisms['stdNa'].plot_kinetics(ax=ax, linestyle='--')

Got data for v in range -100.0 to 100.0 at 37°C
Got data for v in range -100.0 to 100.0 at 23°C

WARNING: Mechanism Na is not inserted in any domain and will be removed.
(model.py, line 637)

Saved content to examples/Toy/biophys/mod/stdNa.mod
Got data for v in range -100.0 to 100.0 at 37°C

Let’s do the same for the potassium channel:

fig, ax = plt.subplots(1, 2, figsize=(10, 5))
model.mechanisms['Kv'].plot_kinetics(ax=ax)
model.standardize_channel('Kv');
model.mechanisms['stdKv'].plot_kinetics(ax=ax, linestyle='--')

Got data for v in range -100.0 to 100.0 at 37°C
Got data for v in range -100.0 to 100.0 at 23°C

WARNING: Mechanism Kv is not inserted in any domain and will be removed.
(model.py, line 637)

Saved content to examples/Toy/biophys/mod/stdKv.mod
Got data for v in range -100.0 to 100.0 at 37°C

We can plot the voltage trace of the standardized model to compare it with the trace produced by the original model.

model.run(300)

fig, ax = plt.subplots(1, 1, figsize=(10, 4))
ax.plot(t, v_orig, color='C0', linestyle='-', zorder=0, lw=1, label='original')
model.simulator.plot_voltage(ax=ax, linestyle='--', color='C1')