Morphology reduction (Notebook)

In this tutorial, we will learn how to reduce a subtree of a neuron to a single equivalent cylinder. We will also see how to map the parameters of the original segments to the reduced segments.

We will follow the analytical impedance-based approach proposed by Amsalem et al., 2020. We have integrated this reduction algorithm in DendroTweaks and extended its functionality to allow for a continuum of morphology reduction levels, bridging detailed and “ball-and-stick”-like models. You can select any section of the cell and reduce its subtree, which allows for any intermediate level of detail to be achieved.

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

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

Loading a model

Let’s create a model:

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

Path to model: examples/Park_2019

model = dd.Model(path_to_model)

model.list_morphologies()

['Park_2019_sorted', 'Park_2019']

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).

First, we need to visualize the section tree of our model.

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

Next we add biophysical mechanisms to the model.

model.list_biophys()

['Park_2019_full']

model.load_biophys('Park_2019_full')

Saved content to examples/Park_2019/biophys/python/Km.py
Saved content to examples/Park_2019/biophys/python/Na.py
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/Kv.py
Saved content to examples/Park_2019/biophys/python/Ka.py

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}

We will visualize and compare the parameter distributions before and after model reduction. Let’s see how the distributions look before reduction:

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

Let’s also run a simulation before reduction so that we can compare the voltage traces berfore and after:

model.list_stimuli()

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

soma = model.sec_tree.soma
model.load_stimuli('depolarizing_current')

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

model.run(1000)

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)]

Morphology reduction

Now, let’s consider a scenario where we want to reduce the entire apical subtree of the cell to a single equivalent cylinder. To accomplish this reduction, we need to select a section that serves as the root of the subtree we want to merge. As you might notice at the section tree plot above, section 9 is the root of the entire apical dendrite. Let’s select this section for reduction:

root = model.sec_tree[9]

We can visualize the geometry of the section before reduction:

fig, ax = plt.subplots(figsize=(10, 2))
root.plot_radii(ax=ax)

<Axes: title={'center': 'Radius Distribution - Section 9 (apic)'}, xlabel='Distance (µm)', ylabel='Radius (µm)'>

Now we can use the reduce_subtree method to run the reduction algorithm, calculate the properties of the equivalent cylinder, and replace the entire subtree with a single section. We set the fit parameter to True to fit a polynomial to the channel distributions. This allows us to represent the channel distribution efficiently with just a few coefficients instead of storing values for each segment explicitly.

Note that the algorithm doesn’t currently support reduction with stimuli or recordings in the subtree. In our case, this isn’t an issue since we record and stimulate only at the soma.

data = model.reduce_subtree(root, fit=True)

CableParams(length=764.4908970853801, diam=2.5167725013593345, space_const=831.9329527514925, cm=2.0, rm=11000.0, ra=100.0, e_pas=-79.0, electrotonic_length=0.9189333041285863)
Sorted SectionTree(root=NeuronSection(idx=0), num_nodes=45).
Sorted PointTree(root=Point(idx=0), num_nodes=2048).
Sorted SegmentTree(root=NeuronSegment(idx=0), num_nodes=219).
Sorted SegmentTree(root=NeuronSegment(idx=0), num_nodes=62).
Sorted PointTree(root=Point(idx=0), num_nodes=396).
Interpolated for ids []

WARNING: Group distal_apical is empty and will be removed.
(model.py, line 776)
WARNING: Fitting failed for parameter gbar_CaLVA with the provided tolerance.
Using the last valid fit (degree=6). Maximum residual: 1.7407861052657653e-07
(model.py, line 1511)
WARNING: Fitting failed for parameter gbar_Ka with the provided tolerance.
Using the last valid fit (degree=6). Maximum residual: 0.002857036610448546
(model.py, line 1511)
WARNING: Domain apic is empty and will be removed.
(model.py, line 751)
WARNING: Group apical is empty and will be removed.
(model.py, line 776)
WARNING: Group proximal_apical is empty and will be removed.
(model.py, line 776)

Let’s visualize the mapping from original to reduced segments. Each point on the plot represents an original segment, with the x-axis showing its normalized distance along the subtree. The colors indicate how the original segments map to the segment of the reduced section.

Reduced model

segs_to_reduced_segs = data['segs_to_reduced_segs']
segs_to_locs = data['segs_to_locs']
fig, ax = plt.subplots(figsize=(10, 2))
reduced_segs = root.segments
reduced_seg_to_color = {seg: plt.cm.tab10.colors[i] for i, seg in enumerate(reduced_segs)}
ax.set_ylim(-1, 1)
for seg, loc in segs_to_locs.items():
    plt.scatter(loc, [0], s=10, color=reduced_seg_to_color[segs_to_reduced_segs[seg]])

norm_centers = np.round(np.array(root.seg_centers)/root.L, 3)
ax.set_xticks(norm_centers);
ax.set_xticklabels(norm_centers);
ax.set_yticks([]);

We have 9 segments in the reduced section. Let’s now plot the section’s geometry as we did before:

fig, ax = plt.subplots(figsize=(10, 2))
root.plot_radii(ax=ax)

<Axes: title={'center': 'Radius Distribution - Section 9 (reduced_0)'}, xlabel='Distance (µm)', ylabel='Radius (µm)'>

Let’s examine how the domains have changed. You’ll notice a new reduced_0 domain containing only the reduced section. There’s also a new group with the same name that includes only the segments of the reduced section. This allows the section to have unique distribution parameters that represent the averaged distribution within the original subtree.

model.domains_to_mechs

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

model.groups

{'all': SegmentGroup("all", domains=['axon', 'dend', 'soma', 'reduced_0']),
 'axonal': SegmentGroup("axonal", domains=['axon']),
 'dendritic': SegmentGroup("dendritic", domains=['dend']),
 'somatic': SegmentGroup("somatic", domains=['soma']),
 'dendritic_thin': SegmentGroup("dendritic_thin", domains=['dend'], section_diam(None, 0.8)),
 'dendritic_thick': SegmentGroup("dendritic_thick", domains=['dend'], section_diam(0.8, None)),
 'basal': SegmentGroup("basal", domains=['dend']),
 'reduced_0': SegmentGroup("reduced_0", domains=['reduced_0'])}

Let’s examine how the model parameters have changed. Notice that we have a new reduced_0 group and its associated polynomial distributions with their coefficients.

model.params

{'cm': {'all': constant({'value': 2}),
  'somatic': constant({'value': 1}),
  'reduced_0': constant({'value': 2.0})},
 'Ra': {'all': constant({'value': 100}),
  'reduced_0': constant({'value': 100.0})},
 'gbar_Leak': {'all': constant({'value': 9.09090909090909e-05}),
  'reduced_0': constant({'value': 9.090909090909092e-05})},
 'e_Leak': {'all': constant({'value': -79}),
  'reduced_0': constant({'value': -79.0})},
 'depth_CaDyn': {'all': constant({'value': 0.1}),
  'reduced_0': constant({'value': 0.10000000000000003})},
 'taur_CaDyn': {'all': constant({'value': 50}),
  'reduced_0': constant({'value': 50.0})},
 'cainf_CaDyn': {'all': constant({'value': 0.0001}),
  'reduced_0': constant({'value': 0.00010000000000000003})},
 'gamma_CaDyn': {'all': constant({'value': 1}),
  'reduced_0': constant({'value': 1.0})},
 'kt_CaDyn': {'all': constant({'value': 0}),
  'reduced_0': constant({'value': -0.0})},
 'kd_CaDyn': {'all': constant({'value': 0}),
  'reduced_0': constant({'value': -0.0})},
 'gbar_CaHVA': {'all': constant({'value': 0.0}),
  'somatic': constant({'value': 5e-06}),
  'basal': linear({'slope': 1e-08, 'intercept': 5e-06}),
  'reduced_0': polynomial({'coeffs': array([-8.26293184e-21,  1.97195076e-17, -1.82485957e-14,  8.19306560e-12,
         -1.79431090e-09,  1.57125582e-07,  7.30059330e-07])})},
 'eca': {'all': constant({'value': 140}),
  'reduced_0': constant({'value': 113.2472408069})},
 'gbar_CaLVA': {'all': constant({'value': 0.0}),
  'somatic': constant({'value': 3e-05}),
  'basal': linear({'slope': 6e-08, 'intercept': 3e-05}),
  'reduced_0': polynomial({'coeffs': array([-4.95775910e-20,  1.18317045e-16, -1.09491574e-13,  4.91583936e-11,
         -1.07658654e-08,  9.42753492e-07,  4.38035598e-06])})},
 'gbar_KCa': {'all': constant({'value': 0.0}),
  'somatic': constant({'value': 0.00021}),
  'dendritic': constant({'value': 0.00021}),
  'reduced_0': constant({'value': 0.0002100000000000001})},
 'ek': {'all': constant({'value': -80}),
  'reduced_0': constant({'value': -80.0})},
 'gbar_Ka': {'all': constant({'value': 0.0}),
  'somatic': constant({'value': 0.0054}),
  'dendritic_thin': constant({'value': 0.108}),
  'dendritic_thick': constant({'value': 0.0108}),
  'reduced_0': polynomial({'coeffs': array([-1.50642424e-17,  2.94201489e-14, -1.94016044e-11,  4.38844154e-09,
         -2.29108662e-08,  3.91907843e-05,  7.38410860e-02])})},
 'gbar_Km': {'all': constant({'value': 0.0}),
  'somatic': constant({'value': 0.0002794}),
  'dendritic': constant({'value': 0.000127}),
  'reduced_0': constant({'value': 0.00012700000000000008})},
 'v12_Km': {'all': constant({'value': -30}),
  'reduced_0': constant({'value': -30.0})},
 'q_Km': {'all': constant({'value': 9}),
  'reduced_0': constant({'value': 9.0})},
 'gbar_Kv': {'all': constant({'value': 0.0}),
  'somatic': constant({'value': 0.005}),
  'dendritic': constant({'value': 0.00015}),
  'reduced_0': constant({'value': 0.00015000000000000007})},
 'v12_Kv': {'all': constant({'value': 25}),
  'reduced_0': constant({'value': 25.0})},
 'q_Kv': {'all': constant({'value': 9}),
  'reduced_0': constant({'value': 9.0})},
 'gbar_Na': {'all': constant({'value': 0.0}),
  'somatic': constant({'value': 0.0505}),
  'dendritic': constant({'value': 0.0303}),
  'reduced_0': constant({'value': 0.030300000000000007})},
 'Rma_Na': {'all': constant({'value': 0.182}),
  'reduced_0': constant({'value': 0.18200000000000008})},
 'Rmb_Na': {'all': constant({'value': 0.14}),
  'reduced_0': constant({'value': 0.14000000000000007})},
 'v12m_Na': {'all': constant({'value': -30}),
  'reduced_0': constant({'value': -30.0})},
 'qm_Na': {'all': constant({'value': 9.8}),
  'reduced_0': constant({'value': 9.8})},
 'Rhb_Na': {'all': constant({'value': 0.0091}),
  'reduced_0': constant({'value': 0.009100000000000004})},
 'Rha_Na': {'all': constant({'value': 0.024}),
  'reduced_0': constant({'value': 0.024000000000000007})},
 'v12ha_Na': {'all': constant({'value': -45}),
  'reduced_0': constant({'value': -45.0})},
 'v12hb_Na': {'all': constant({'value': -70}),
  'reduced_0': constant({'value': -70.0})},
 'qh_Na': {'all': constant({'value': 5}),
  'reduced_0': constant({'value': 5.0})},
 'v12hinf_Na': {'all': constant({'value': -60}),
  'reduced_0': constant({'value': -60.0})},
 'qhinf_Na': {'all': constant({'value': 6.2}),
  'reduced_0': constant({'value': 6.2})},
 'ena': {'all': constant({'value': 60}),
  'reduced_0': constant({'value': 60.0})}}

Let’s plot the parameter distributions in the reduced model:

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

Finally, we can plot the section tree to see how the cell looks after reduction

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

Let’s run a simulation using the reduced model and compare the results to the original model output.

Comparing simulation results

model.run(1000)

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')

Let’s slightly increase the current amplitude and run the simulation again to match the original model:

print(f"Initial amplitude: {model.iclamps[soma(0.5)].amp} nA")
model.iclamps[soma(0.5)].amp= 0.172
print(f"New amplitude: {model.iclamps[soma(0.5)].amp} nA")

Initial amplitude: 0.162 nA
New amplitude: 0.172 nA

model.run(1000)

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')