Quick start

This quick start tutorial will guide you through the basic steps of creating single-cell biophysical neuronal models in DendroTweaks.

We will use the standalone DendroTweaks Python package to build a detailed biophysical model of a L2/3 pyramidal neuron from the mouse visual cortex (Park et al., 2019). We will start by loading and examining morphological reconstruction from an SWC file. Next, we will add membrane mechanisms, such as ion channels, from MOD files and insert them into cell domains (e.g., apical, basal). We will then learn how to adjust different model parameters and run simulations. Finally, we will use built-in validation protocols to examine the responses of our model in order to be able to compare them to experimental data.

Alternatively, you can run Park_2019 model using the GUI provided on our online platform at https://dendrotweaks.dendrites.gr.

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)

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

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

import dendrotweaks as dd
dd.__version__

'0.4.2'

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

Organizing model data

Before diving into model building, we need to establish a well-organized directory structure for our model data. First, we create the examples directory that will store all of our models and download the example models and their corresponding data from the Github repository.

/path/to/examples/
    ├── Park_2019/
    ├── Poirazi_2003/
    ...

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

Downloading example data...
Downloading examples from https://github.com/Poirazi-Lab/DendroTweaks/archive/refs/heads/main.zip...
Extracting examples to examples
Examples downloaded successfully to examples/.

os.listdir('examples/')

['Smith_2013', 'Park_2019', 'Hay_2011', 'Test', 'Poirazi_2003']

To create a model, we must specify the path to the folder containing the model files.

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

Path to model: examples/Park_2019

Each model requires two essential components: neuronal morphology (stored as SWC files in the morphology folder) and membrane mechanisms (ion channels defined as MOD files in the biophys/mod folder).

/path/to/examples/
    ...
    └── Park_2019/
        ├── morphology/
        |  └── Park_2019.swc
        └── biophys/
           └── mod/
              ├── Nav.mod
              ...
              └── Kv.mod

Tip

Don’t worry about compiling the MOD files manually - DendroTweaks will handle this step automatically later in the process.

Creating a model

With our directory structure in place, we can now create a model object. This object serves as the central interface for managing all model components, including morphology, kinetics and distribution of membrane mechanisms, and stimuli.

model = dd.Model(path_to_model)

When we instantiate the model by specifying the path to our model directory, DendroTweaks automatically sets up additional required folders. It creates a Default folder containing standard mechanisms and a Templates folder with JINJA templates (e.g. for converting MOD files to Python classes). These folders are created only once and are shared across all models.

Within your model folder, it generates a stimuli directory for stimulation protocols and a biophys/python folder for storing Python files generated from MOD files.

This structure provides the foundation for model development. As we progress, we will add more files to these directories, but the basic framework is now ready. You can learn more about the directory structure in the tutorial on loading and saving models.

Morphology

We will proceed by loading the morphology of the cell from an SWC file. First, we list the available morphologies in the morphology subfolder.

model.list_morphologies()

['Park_2019_sorted', 'Park_2019']

We can load a specific morphology using the load_morphology method.

model.load_morphology(file_name='Park_2019')

Sorted PointTree(root=Point(idx=0), num_nodes=2214).
Extended 44 nodes.
Sorted PointTree(root=Point(idx=0), num_nodes=2258).

In DendroTweaks, we represent neuronal morphologies as hierarchical tree graphs. A single neuron can be represented in three distinct ways: the point tree, the section tree, and the segment tree. Each tree graph contains nodes that correspond to different elements: points from a morphology reconstruction (SWC file), sections, or spatial segments. The edges between nodes show parent-child relationships. These three representations capture different levels of abstraction: geometry (point tree), topology (section tree), and spatial discretization (segment tree).

These tree-based representations work together to provide complementary views of the cell’s structural organization and spatial relationships. They are interconnected, with each section linking to both its geometry (points) and its segmentation (segments), and vice versa.

For more information about representing neuronal morphology with tree graphs, see the corresponding tutorial.

model.point_tree, model.sec_tree, model.seg_tree

(PointTree(root=Point(idx=0), num_nodes=2258),
 SectionTree(root=NeuronSection(idx=0), num_nodes=52),
 SegmentTree(root=NeuronSegment(idx=0), num_nodes=124))

We can visualize the section tree using the plot method:

fig, ax = plt.subplots(figsize=(10, 10))
model.sec_tree.plot(ax, 
                    show_points=True,
                    show_lines=True,
                    annotate=True)

We can calculate some important morphometric statistics for distinct anatomical and functional domains of the cell with the following function from the analysis subpackage:

from dendrotweaks.analysis import calculate_domain_statistics

calculate_domain_statistics(model, domain_names=['apic', 'soma'])

{'apic': {'N_sections': 43,
  'N_bifurcations': 21,
  'N_terminations': 22,
  'depth': {'min': 1,
   'max': 8,
   'counts': {1: 1, 2: 2, 3: 4, 4: 8, 5: 8, 6: 6, 7: 8, 8: 6}},
  'diam': {'min': 0.38, 'max': 2.06, 'mean': 0.82, 'std': 0.45},
  'length': {'min': 4.51, 'max': 154.69, 'mean': 58.31, 'std': 35.75},
  'area': {'min': 17.6, 'max': 271.28, 'mean': 130.52, 'std': 74.2},
  'total_length': 2507.17,
  'total_area': 5612.19},
 'soma': {'N_sections': 1,
  'N_bifurcations': 0,
  'N_terminations': 0,
  'depth': {'min': 0, 'max': 0, 'counts': {0: 1}},
  'diam': {'min': 20.44, 'max': 20.44, 'mean': 20.44, 'std': nan},
  'length': {'min': 20.44, 'max': 20.44, 'mean': 20.44, 'std': nan},
  'area': {'min': 1312.27, 'max': 1312.27, 'mean': 1312.27, 'std': nan},
  'total_length': 20.44,
  'total_area': 1312.27}}

Segmentation

Computer simulations always involve approximating a continuous system as one that is discrete in space and time. Thus, we need to set the spatial discretization of the model. The spatial discretization determines the number of segments in a section.

The optimal number of segments can be calculated based on the frequency-dependent length constant and the spatial discretization coefficient (d_lambda). To learn the underlying equations see the corresponding tutorial. Spatial discretization also depends on the membrane capacitance and axial resistance of the model. The seg_tree we created before was created using default NEURON’s values for capacitance and membrane resistance. In our model, we want to change those values to be more biologically relevant for our neuron.

model.set_param('cm', value=2)
model.set_param('cm', group_name='somatic', value=1)
model.set_param('Ra', value=100)

model.params

{'cm': {'all': constant({'value': 2}), 'somatic': constant({'value': 1})},
 'Ra': {'all': constant({'value': 100})}}

We should interpret the output as follows: the specific membrane capacitance cm is set to a constant value of 2 \(\mu F/cm^2\) for all segments of the cell apart from the soma, where it is set to 1 \(\mu F/cm^2\). The axial resistance Ra is set to 100 \(Ohm \cdot cm\) for all segments of the cell.

We will explore how to set the model parameters in more detail later in the tutorial. For now, we have done just enough to update the spatial discretization.

print(f'Number of segments using the default parameters: {len(model.seg_tree)}')
model.set_segmentation(d_lambda=0.1)
print(f'Number of segments using the new parameters: {len(model.seg_tree)}')

Number of segments using the default parameters: 124
Number of segments using the new parameters: 246

We can visualize a section of the section tree using the plot method. This displays the section’s geometry, with bars in the bottom-right plot corresponding to the segments.

sec = model.sec_tree[22]
fig, ax = plt.subplots(2, 2, figsize=(12, 6))
sec.plot(ax=ax, plot_parent=True)
plt.tight_layout()

Membrane mechanisms

After defining the neuronal morphology, we now need to specify biophysical properties. The biophysical properties of our model depend on the mechanisms present in the membrane. These mechanisms are defined in MOD files as sets of equations and parameters, which are compiled and loaded into NEURON.

Adding mechanisms

First, let’s add the default mechanisms to the model, which include the leak channel, calcium dynamics and synaptic mechanisms:

model.add_default_mechanisms(recompile=True)

However, for most of the models, we need to add user-defined mechanisms. We can create a mechanism object from a MOD file using the add_mechanism method. We can list the available mechanisms in the biophys/mod directory:

model.list_mechanisms()

['KCa', 'Ka', 'Km', 'Na', 'CaLVA', 'CaHVA', 'Kv']

We can add these user-defined mechanisms to the model:

for mech_name in ['CaHVA', 'CaLVA', 'KCa', 'Ka', 'Km', 'Kv', 'Na']:
    model.add_mechanism(mech_name, recompile=True)

Saved content to examples/Park_2019/biophys/python/CaHVA.py
Saved content to examples/Park_2019/biophys/python/CaLVA.py
Saved content to examples/Park_2019/biophys/python/KCa.py
Saved content to examples/Park_2019/biophys/python/Ka.py
Saved content to examples/Park_2019/biophys/python/Km.py
Saved content to examples/Park_2019/biophys/python/Kv.py
Saved content to examples/Park_2019/biophys/python/Na.py

With these commands, we create Python objects from MOD files, that we will use to visualize chanel kinetics. We also compile and load the MOD files, making them available in NEURON.

To see all mechanisms available in the model, we can use the mechanisms attribute.

model.mechanisms

{'Leak': <Mechanism(Leak)>,
 'CaDyn': <Mechanism(CaDyn)>,
 'CaHVA': <Mechanism(CaHVA)>,
 'CaLVA': <Mechanism(CaLVA)>,
 'KCa': <Mechanism(KCa)>,
 'Ka': <Mechanism(Ka)>,
 'Km': <Mechanism(Km)>,
 'Kv': <Mechanism(Kv)>,
 'Na': <Mechanism(Na)>}

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}

Warning

Note that the parameters stored within the mechanisms are the default values from the MOD files. The actual values of the parameters used for the simulation are stored in the model object!

We can view the global parameters of the model with the params attribute:

model.params

{'cm': {'all': constant({'value': 2}), 'somatic': constant({'value': 1})},
 'Ra': {'all': constant({'value': 100})}}

Warning

Note that so far we have only loaded the mechanisms without actually inserting them into the membrane. Therefore, the parameters of these mechanisms are not yet included in the model.params dictionary. In the next step, we will insert the mechanisms into the membrane.

Inserting mechanisms to specific domains

In DendroTweaks, membrane mechanisms are mapped to morphological domains. A domain is a region of a neuron distinguished by its anatomical or functional properties. In a typical pyramidal cell model we have the following domains: soma, axon, basal dendrites, apical dendrite (further subdivided into trunk, tuft, and oblique dendrites).

In DendroTweaks, a domain represents a collection of sections that share the same properties. We can view the domains of the model with the domains attribute.

model.domains

{'apic': <Domain(apic, 43 sections)>,
 'axon': <Domain(axon, 1 sections)>,
 'dend': <Domain(dend, 7 sections)>,
 'soma': <Domain(soma, 1 sections)>}

Following the mechanism upload in the previous step, we need to insert these mechanisms into specific domains. For this example, we insert all available mechanisms into every domain except the axon, which receives only the leak channel.

for domain_name in model.domains:
    if domain_name == 'axon':
        model.insert_mechanism('Leak', 'axon', distribute=False)
        continue
    for mech_name in model.mechanisms:
        model.insert_mechanism(mech_name, domain_name, distribute=False)

We can verify which mechanisms are inserted in each domain with the domains_to_mechs attribute.

model.domains_to_mechs

{'apic': {'CaDyn', 'CaHVA', 'CaLVA', 'KCa', 'Ka', 'Km', 'Kv', 'Leak', 'Na'},
 'axon': {'Leak'},
 'dend': {'CaDyn', 'CaHVA', 'CaLVA', 'KCa', 'Ka', 'Km', 'Kv', 'Leak', 'Na'},
 'soma': {'CaDyn', 'CaHVA', 'CaLVA', 'KCa', 'Ka', 'Km', 'Kv', 'Leak', 'Na'}}

To examine the parameters of the inserted mechanisms, we can use the mechs_to_params attribute.

model.mechs_to_params

{'Independent': ['cm', 'Ra', 'eca', 'ek', 'ena'],
 'Leak': ['gbar_Leak', 'e_Leak'],
 'CaDyn': ['depth_CaDyn',
  'taur_CaDyn',
  'cainf_CaDyn',
  'gamma_CaDyn',
  'kt_CaDyn',
  'kd_CaDyn'],
 'CaHVA': ['gbar_CaHVA'],
 'CaLVA': ['gbar_CaLVA'],
 'KCa': ['gbar_KCa'],
 'Ka': ['gbar_Ka'],
 'Km': ['gbar_Km', 'v12_Km', 'q_Km'],
 'Kv': ['gbar_Kv', 'v12_Kv', 'q_Kv'],
 'Na': ['gbar_Na',
  'Rma_Na',
  'Rmb_Na',
  'v12m_Na',
  'qm_Na',
  'Rhb_Na',
  'Rha_Na',
  'v12ha_Na',
  'v12hb_Na',
  'qh_Na',
  'v12hinf_Na',
  'qhinf_Na']}

Some parameters, such as specific membrane capacitance cm and axial resistance Ra, do not belong to any mechanism. These independent parameters are grouped under an “Independent” pseudo-mechanism for interface consistency. These parameters are available in each domain by default.

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 model using the params attribute.

model.df_params

	Mechanism	Parameter	Group	Distribution	Distribution params
0	Independent	cm	all	constant	{'value': 2}
1	Independent	cm	somatic	constant	{'value': 1}
2	Independent	Ra	all	constant	{'value': 100}
3	Independent	eca	all	constant	{'value': 140}
4	Independent	ek	all	constant	{'value': -77}
5	Independent	ena	all	constant	{'value': 50}
6	Leak	gbar_Leak	all	constant	{'value': 0.0}
7	Leak	e_Leak	all	constant	{'value': -70}
8	CaDyn	depth_CaDyn	all	constant	{'value': 0.1}
9	CaDyn	taur_CaDyn	all	constant	{'value': 80}
10	CaDyn	cainf_CaDyn	all	constant	{'value': 0.0001}
11	CaDyn	gamma_CaDyn	all	constant	{'value': 0.05}
12	CaDyn	kt_CaDyn	all	constant	{'value': 0.0}
13	CaDyn	kd_CaDyn	all	constant	{'value': 0.0}
14	CaHVA	gbar_CaHVA	all	constant	{'value': 0.0}
15	CaLVA	gbar_CaLVA	all	constant	{'value': 0.0}
16	KCa	gbar_KCa	all	constant	{'value': 0.0}
17	Ka	gbar_Ka	all	constant	{'value': 0.0}
18	Km	gbar_Km	all	constant	{'value': 0.0}
19	Km	v12_Km	all	constant	{'value': -30}
20	Km	q_Km	all	constant	{'value': 9}
21	Kv	gbar_Kv	all	constant	{'value': 0.0}
22	Kv	v12_Kv	all	constant	{'value': 25}
23	Kv	q_Kv	all	constant	{'value': 9}
24	Na	gbar_Na	all	constant	{'value': 0.0}
25	Na	Rma_Na	all	constant	{'value': 0.182}
26	Na	Rmb_Na	all	constant	{'value': 0.14}
27	Na	v12m_Na	all	constant	{'value': -30}
28	Na	qm_Na	all	constant	{'value': 9.8}
29	Na	Rhb_Na	all	constant	{'value': 0.0091}
30	Na	Rha_Na	all	constant	{'value': 0.024}
31	Na	v12ha_Na	all	constant	{'value': -45}
32	Na	v12hb_Na	all	constant	{'value': -70}
33	Na	qh_Na	all	constant	{'value': 5}
34	Na	v12hinf_Na	all	constant	{'value': -60}
35	Na	qhinf_Na	all	constant	{'value': 6.2}

As we can see, all parameters are set to their default values (from the MOD files) across all segments of the cell, with some parameters initialized to 0.0. Before running the simulation, we need to set these parameters to more realistic values, which we will learn how to do in the next step.

Model parameters

In real neurons, properties such as ion channel density vary across different regions of the cell. To distribute parameters across the cell, we need to specify where and how the parameter will be distributed.

To select the segments where a given distribution will be applied, we use segment groups. A segment group is a collection of segments that meet certain criteria, such as the diameter or distance from the soma.

To define how the parameter will be distributed, we use distribution functions. A distribution function takes a segment’s distance from the soma as input and returns the parameter value at that distance. The figure below shows an example of a segment group for the apical nexus region and a Gaussian distribution function for a parameter, such as ion channel conductance.

Defining segment groups

By default DendroTweaks creates a group for each domain and the all group for the entire cell.

model.groups

{'all': SegmentGroup("all", domains=['apic', 'axon', 'dend', 'soma']),
 'apical': SegmentGroup("apical", domains=['apic']),
 'axonal': SegmentGroup("axonal", domains=['axon']),
 'dendritic': SegmentGroup("dendritic", domains=['dend']),
 'somatic': SegmentGroup("somatic", domains=['soma'])}

We will need more specific groups to reproduce biologicaly realistic distribution of ion channels in the cell. Available values for the select_by argument are:

• diam - diameter of the segment

• section_diam - diameter at the center of the section to which the segment belongs

• distance - distance of the segment center from the soma center

• domain_distance - distance of the segment center to the closest parent segment in a different domain

model.add_group(
    'dendritic_thin', domains=['dend', 'apic'],
    select_by='section_diam', max_value=0.8
)
model.add_group(
    'dendritic_thick', domains=['dend', 'apic'],
    select_by='section_diam', min_value=0.8
)

model.add_group(
    'proximal_apical', domains=['apic'],
    select_by='distance', max_value=260
)
model.add_group(
    'distal_apical', domains=['apic'],
    select_by='distance', min_value=260
)

model.add_group('basal', domains=['dend'])

model.add_group('dendritic', domains=['dend', 'apic'])

Distributing ion channels

To define how the parameter will be distributed, we use distribution functions. A distribution function takes a segment’s distance from the soma as input and returns the parameter value at that distance.

We can set the values of the parameters for the mechanisms inserted in the model using the set_param method, specifying the group name and the distribution type.

# Leak
Rm = 11000
model.set_param('gbar_Leak', value=1/Rm)

# Na
model.set_param('gbar_Na', group_name='somatic', value=0.0505)
model.set_param('gbar_Na', group_name='dendritic', value=0.0303)

# # Kv
model.set_param('gbar_Kv', group_name='somatic', value=0.005)
model.set_param('gbar_Kv', group_name='dendritic', value=0.000_15)

# # Km
model.set_param('gbar_Km', group_name='somatic', value=0.0002794)
model.set_param('gbar_Km', group_name='dendritic', value=0.000127)

# Ka
model.set_param('gbar_Ka', group_name='somatic', value=0.0054)
model.set_param('gbar_Ka', group_name='dendritic_thin', value=0.108)
model.set_param('gbar_Ka', group_name='dendritic_thick', value=0.0108)

# Kca
model.set_param('gbar_KCa', group_name='somatic', value=0.000_21)
model.set_param('gbar_KCa', group_name='dendritic', value=0.000_21)

The above channels have simple constant distribution with some difference between groups. Note that we will use more complex distributions for the calcium channels:

#CaHVA
model.set_param('gbar_CaHVA', group_name='somatic', value=0.000_005)
model.set_param('gbar_CaHVA', group_name='basal', 
    distr_type='linear', slope=0.000_000_01, intercept=0.000_005)
model.set_param('gbar_CaHVA', group_name='proximal_apical', 
    distr_type='sinusoidal', amplitude=0.000_004923, frequency=0.008758, phase=0.8656)
model.set_param('gbar_CaHVA', group_name='distal_apical', value=0.000_002)

# CaLVA
model.set_param('gbar_CaLVA', group_name='somatic', value=0.000_03)
model.set_param('gbar_CaLVA', group_name='basal', 
    distr_type='linear', slope=0.000_000_06, intercept=0.000_03)
model.set_param('gbar_CaLVA', group_name='proximal_apical', 
    distr_type='sinusoidal', amplitude=0.000_029538, frequency=0.008758, phase=0.8656)
model.set_param('gbar_CaLVA', group_name='distal_apical', value=0.000_012)

The following distribution functions are available, along with their expected parameters:

• constant: Requires a value parameter.

• linear: Requires slope and intercept parameters.

• exponential: Requires vertical_shift, scale_factor, growth_rate, and horizontal_shift parameters.

• sigmoid: Requires vertical_shift, scale_factor, growth_rate, and horizontal_shift parameters.

• sinusoidal: Requires amplitude, frequency, and phase parameters.

• gaussian: Requires amplitude, mean, and std parameters.

We can also set up the parameters for internal calcium dynamics:

# Internal calcium dynamics
model.set_param('depth_CaDyn', value=0.1)
model.set_param('taur_CaDyn', value=50)
model.set_param('cainf_CaDyn', value=1e-4)
model.set_param('gamma_CaDyn', value=1)
model.set_param('kt_CaDyn', value=0)
model.set_param('kd_CaDyn', value=0)

We can visualize the distribution of ion channels conductances in our model as a function of distance from the soma.

for param_name in model.conductances:
    fig, ax = plt.subplots(figsize=(10, 2))
    model.plot_param(param_name, ax=ax)

We can also set other parameters, such as reversal potentials, temperature, and initial membrane potential.

model.set_param('e_Leak', value=-79)
if 'ena' in model.params:
    model.set_param('ena', value=60)
if 'ek' in model.params:
    model.set_param('ek', value=-80)
if 'eca' in model.params:
    model.set_param('eca', value=140)
model.set_param('temperature', value=37)
model.set_param('v_init', value=-79)

# model.distribute_all()

Now, we can access the model parameters again and see that the values have been updated.

model.df_params

	Mechanism	Parameter	Group	Distribution	Distribution params
0	Independent	cm	all	constant	{'value': 2}
1	Independent	cm	somatic	constant	{'value': 1}
2	Independent	Ra	all	constant	{'value': 100}
3	Independent	eca	all	constant	{'value': 140}
4	Independent	ek	all	constant	{'value': -80}
5	Independent	ena	all	constant	{'value': 60}
6	Leak	gbar_Leak	all	constant	{'value': 9.09090909090909e-05}
7	Leak	e_Leak	all	constant	{'value': -79}
8	CaDyn	depth_CaDyn	all	constant	{'value': 0.1}
9	CaDyn	taur_CaDyn	all	constant	{'value': 50}
10	CaDyn	cainf_CaDyn	all	constant	{'value': 0.0001}
11	CaDyn	gamma_CaDyn	all	constant	{'value': 1}
12	CaDyn	kt_CaDyn	all	constant	{'value': 0}
13	CaDyn	kd_CaDyn	all	constant	{'value': 0}
14	CaHVA	gbar_CaHVA	all	constant	{'value': 0.0}
15	CaHVA	gbar_CaHVA	somatic	constant	{'value': 5e-06}
16	CaHVA	gbar_CaHVA	basal	linear	{'slope': 1e-08, 'intercept': 5e-06}
17	CaHVA	gbar_CaHVA	proximal_apical	sinusoidal	{'amplitude': 4.923e-06, 'frequency': 0.008758...
18	CaHVA	gbar_CaHVA	distal_apical	constant	{'value': 2e-06}
19	CaLVA	gbar_CaLVA	all	constant	{'value': 0.0}
20	CaLVA	gbar_CaLVA	somatic	constant	{'value': 3e-05}
21	CaLVA	gbar_CaLVA	basal	linear	{'slope': 6e-08, 'intercept': 3e-05}
22	CaLVA	gbar_CaLVA	proximal_apical	sinusoidal	{'amplitude': 2.9538e-05, 'frequency': 0.00875...
23	CaLVA	gbar_CaLVA	distal_apical	constant	{'value': 1.2e-05}
24	KCa	gbar_KCa	all	constant	{'value': 0.0}
25	KCa	gbar_KCa	somatic	constant	{'value': 0.00021}
26	KCa	gbar_KCa	dendritic	constant	{'value': 0.00021}
27	Ka	gbar_Ka	all	constant	{'value': 0.0}
28	Ka	gbar_Ka	somatic	constant	{'value': 0.0054}
29	Ka	gbar_Ka	dendritic_thin	constant	{'value': 0.108}
30	Ka	gbar_Ka	dendritic_thick	constant	{'value': 0.0108}
31	Km	gbar_Km	all	constant	{'value': 0.0}
32	Km	gbar_Km	somatic	constant	{'value': 0.0002794}
33	Km	gbar_Km	dendritic	constant	{'value': 0.000127}
34	Km	v12_Km	all	constant	{'value': -30}
35	Km	q_Km	all	constant	{'value': 9}
36	Kv	gbar_Kv	all	constant	{'value': 0.0}
37	Kv	gbar_Kv	somatic	constant	{'value': 0.005}
38	Kv	gbar_Kv	dendritic	constant	{'value': 0.00015}
39	Kv	v12_Kv	all	constant	{'value': 25}
40	Kv	q_Kv	all	constant	{'value': 9}
41	Na	gbar_Na	all	constant	{'value': 0.0}
42	Na	gbar_Na	somatic	constant	{'value': 0.0505}
43	Na	gbar_Na	dendritic	constant	{'value': 0.0303}
44	Na	Rma_Na	all	constant	{'value': 0.182}
45	Na	Rmb_Na	all	constant	{'value': 0.14}
46	Na	v12m_Na	all	constant	{'value': -30}
47	Na	qm_Na	all	constant	{'value': 9.8}
48	Na	Rhb_Na	all	constant	{'value': 0.0091}
49	Na	Rha_Na	all	constant	{'value': 0.024}
50	Na	v12ha_Na	all	constant	{'value': -45}
51	Na	v12hb_Na	all	constant	{'value': -70}
52	Na	qh_Na	all	constant	{'value': 5}
53	Na	v12hinf_Na	all	constant	{'value': -60}
54	Na	qhinf_Na	all	constant	{'value': 6.2}

To learn more about segment groups and parameter distributions, refer to the tutorial on distributing parameters.

We can export the membrane configuration and store it in a JSON file:

model.export_biophys(file_name='Park_2019_full', indent=4)

Simulation and analysis

We will learn how to simulate neuronal activity by applying current stimuli and adding synapses to the neuron model and recording its response. This process mimics experimental electrophysiology where researchers inject current into neurons to study their firing properties.

We will increase the numerical simulation time step from \(0.025 \ ms\) to \(0.01\ ms\), which will speed up the simulations while only slightly reducing precision.

model.simulator.dt = 0.01

Somatic spikes

Let’s begin by examining somatic action potentials. For this we will inject depolarizing current step at the soma and record the somatic voltage.

Let’s first create a new variable for easy access to the soma

soma = model.sec_tree[0]
soma

NeuronSection(idx=0)

We add a recording point at the center of the soma. The loc parameter specifies the location along the section where the recording will be placed. It is a normalized length, with 0.0 representing the start of the section and 1.0 representing the end.

model.add_recording(soma, loc=0.5)

Recording added to sec NeuronSection(idx=0) at loc 0.5.

Now we will apply a current step stimulus to drive the neuron to fire action potentials. This mimics the experimental technique where constant current is injected into a neuron.

We specify the duration of the stimulus in ms, the delay before the stimulus starts, and the amplitude of the stimulus in nanoamperes.

model.add_iclamp(soma, loc=0.5, amp=0.162, delay=50, dur=900)

IClamp added to sec NeuronSection(idx=0) at loc 0.5.

With our recording and stimulus in place, we can now run a simulation.

model.simulator.run(1000);

This runs the simulation for 1000 milliseconds (1 second). After the simulation completes, the voltage data is stored in simulator.recordings['v'] and the corresponding time points in simulator.t.

We will use built-in DendroTweaks functions to plot the voltage trace and extract spike metrics.

from dendrotweaks.analysis import detect_somatic_spikes, plot_somatic_spikes

fig, ax = plt.subplots(1, 1, figsize=(10, 4))
model.simulator.plot_voltage(ax=ax)
spike_data = detect_somatic_spikes(model)
plot_somatic_spikes(spike_data, ax, show_metrics=True)

Detected 7 spikes
Average spike half-width: 0.96 ms
Average spike amplitude: 79.30 mV
Spike frequency: 7.78 Hz

We can export this specific simulation protocol with all the recordings and stimuli in place:

model.export_stimuli(file_name='depolarizing_current')

model.list_stimuli()

['depolarizing_current',
 'attenuation_soma',
 '50_AMPA_apical',
 'hyperpolarizing_current']

Synapses

Now let’s place some AMPA synapses on the apical tree of our neuron. For this we will create a Population object that represents a population of “virtual” neurons that project to our neuron. To create a population we need to specify the type (AMPA) and the number (50) of synapses to be formed and the segments of the model where they will be randomly distributed. We can then update the input parameters such as rate of incoming inputs and kinetic parameters such as rise and decay time constants.

model.remove_all_stimuli()
model.remove_all_recordings()

model.add_recording(soma, loc=0.5)

segments = model.get_segments(['apical'])
model.add_population(
    segments, 
    N=50, 
    syn_type='AMPA'
)
model.update_population_input_params('AMPA_0', rate=30, end=800, noise=1)
model.update_population_kinetic_params('AMPA_0', tau_rise=0.1, tau_decay=2.5)

Recording added to sec NeuronSection(idx=0) at loc 0.5.
{'rate': 30, 'noise': 1, 'start': 100, 'end': 800, 'weight': 1, 'delay': 0}
{'gmax': 0.001, 'tau_rise': 0.1, 'tau_decay': 2.5, 'e': 0}

fig, ax = plt.subplots(1, 1, figsize=(10, 4))
model.run(1000)
model.simulator.plot_voltage(ax=ax)

We can visualize the spike times of the pre-synaptic “virtual” neurons with the following code:

fig, ax = plt.subplots(1, 1, figsize=(10, 4))
spike_times = model.populations['AMPA']['AMPA_0'].spike_times
for i, (syn, times) in enumerate(spike_times.items()):
    ax.plot(times, np.ones_like(times) * i, '.')

ax.set_xlim(0, 1000)
ax.set_xlabel('Time (ms)')
ax.set_ylabel('Synapse ID')
ax.set_title('AMPA Synaptic Activity')

Text(0.5, 1.0, 'AMPA Synaptic Activity')

To learn more about synapses, refer to the corresponding tutorial.

model.export_stimuli(file_name='50_AMPA_apical')

Voltage attenuation in dendrites

Next, we will explore passive and active properties of dendrites. We will start by evaluating how a stimulus amplitude attenuates as the signal propagates from a dendrite to the soma. For this, we will record the activity from several locations along a path from the tip of a dendrite to the soma.

from dendrotweaks.analysis import calculate_voltage_attenuation, plot_voltage_attenuation

model.remove_all_stimuli()
model.remove_all_recordings()

secs = [model.sec_tree[i] for i in [0, 9, 17, 18, 21, 22]]
for sec in secs:
    model.add_recording(sec, loc=0.5)

model.add_iclamp(secs[-1], loc=0.5, amp=-0.01, delay=100, dur=800)

Recording added to sec NeuronSection(idx=0) at loc 0.5.
Recording added to sec NeuronSection(idx=9) at loc 0.5.
Recording added to sec NeuronSection(idx=17) at loc 0.5.
Recording added to sec NeuronSection(idx=18) at loc 0.5.
Recording added to sec NeuronSection(idx=21) at loc 0.5.
Recording added to sec NeuronSection(idx=22) at loc 0.5.
IClamp added to sec NeuronSection(idx=22) at loc 0.5.

fig, ax = plt.subplots(1, 1, figsize=(10, 4))
model.run(1000)
model.simulator.plot_voltage(ax=ax)

We will use calculate_voltage_attenuation to calculate how the voltage deflection from the resting membrane potential drops as we move further away from the injection site at the dendrite. Note that the rightmost point represents the injection site at the dendrite, which is the maximum distance from the soma.

fig, ax = plt.subplots(1, 1, figsize=(5, 5))
attenuation_data = calculate_voltage_attenuation(model)
plot_voltage_attenuation(attenuation_data, ax)

Stimulating segment: NeuronSegment(idx=103)
Recording segments: [NeuronSegment(idx=0), NeuronSegment(idx=53), NeuronSegment(idx=82), NeuronSegment(idx=85), NeuronSegment(idx=98), NeuronSegment(idx=103)]

model.export_stimuli(file_name='attenuation_soma')

Dendritic nonlinearities

Dendrites are not passive cables; they actively process synaptic inputs through various mechanisms like voltage-gated channels and NMDA receptors. These nonlinearities can significantly impact the neuron’s integrative properties and information processing.

from dendrotweaks.analysis import calculate_dendritic_nonlinearity, plot_dendritic_nonlinearity

To study dendritic nonlinearities, we can gradually increase the synaptic input at a dendritic location and observe the resulting voltage response.

To set up the experiment, we need to place a synapse on a dendritic branch and record the voltage at the same location:

model.remove_all_stimuli()
model.remove_all_recordings()

sec = model.sec_tree[22]
model.add_recording(sec, loc=0.5)
segments = [sec(0.5)]

model.add_recording(sec, loc=0.5)
model.add_population(
    segments, 
    N=1, 
    syn_type='AMPA'
)

Recording added to sec NeuronSection(idx=22) at loc 0.5.
Recording added to sec NeuronSection(idx=22) at loc 0.5.

The calculate_dendritic_nonlinearity function will then gradually increase the synaptic weight (equivalent to increasing the number of synapses, given their synchronous activation) and compare the amplitude of the resulting post-synaptic potentials (PSPs) to the linearly increasing amplitude of a unitary PSP:

fig, ax = plt.subplots(1, 2, figsize=(10, 5))
data = calculate_dendritic_nonlinearity(model, duration=300, max_weight=10, n=10)
plot_dendritic_nonlinearity(data, ax=ax)
ax[0].legend(loc='right')

<matplotlib.legend.Legend at 0x7ad8e778f850>

You can learn more in the tutorial on analyzing simulation results.

---

Congratulations! This concludes our quick start tutorial 🎉🎉🎉

We have covered the fundamental steps to create and simulate a single-cell biophysical neuronal model using DendroTweaks and you are now ready to build more complex models.