Calculating surface area of specific regions in a morphological neuron

Description of problem

Hi everyone,

I have been trying to compute the surface area of specific regions of dendritic trees (and soma) of biophysically detailed neurons in order to then be able to build reduced compartmental models (so neurons would only have 2-3 compartments instead of thousands of them as in biophysically detailed models). Here is the paper which I am trying to apply the approach from. There are two separate aspects which I am not sure about and thought that perhaps someone who has experience with this could provide some insights :). Thank you in advance.

What you have already tried

1. I have an .swc file which has the morphological data and I am loading this to create a (Brian) Morphology object which I can use to calculate the relevant areas. Here is my code for calculating the total area of segments which are labelled as ā€˜somaā€™, ā€˜basalā€™ and ā€˜apicalā€™ (=1, 3 and 4 ā€˜typeā€™ in the .swc file). I am not sure if this is the correct way to calculate the surface areas of these regions in the original morphological neuron, as Iā€™ve never worked with such data before and itā€™s my first attempt. Perhaps I am missing something, but any help or insight would be useful on this.

# Create the morphology
morpho = Morphology.from_swc_file(file_path)

# Initialize surface area counters
soma_area = 0
basal_area = 0
apical_area = 0

# Iterate through the morphology sections to calculate surface areas
def calculate_area(section):
    global soma_area, basal_area, apical_area
    total_area = np.sum(section.area)
    if section.type == 'soma': # Soma
        soma_area += total_area
    elif section.type == 'dend': # Basal dendrite
        basal_area += total_area
    elif section.type == 'apic': # Apical dendrite
        apical_area += total_area

# Use a recursive function to go through all sections
def traverse_sections(section):
    calculate_area(section)
    for child in section.children:
        traverse_sections(child)

# Start traversing from the root section
traverse_sections(morpho)

print(f'Soma surface area: {soma_area}')
print(f'Basal dendrite surface area: {basal_area}')
print(f'Apical dendrite surface area: {apical_area}')

2. This is also a conceptual question. When trying to reduce compartmental models, I am not 100% certain on the terminology used to describe regions such as the ā€˜the soma and first 40um of dendritesā€™ and ā€˜all dendrites laying between 40um and about 240um from the somaā€™, cited from the paper I attached above. Do the 40um/240um distances from the soma mean the Euclidean distances from the somatic compartment? If yes, I wrote a script (below) that tries to calculate the total surface area of all segments within specific boundaries based on the Euclidean distance from the soma. But again, I am not sure if this is the correct way to calculate this or if my approach is fundamentally flawed for what I am trying to accomplish. Perhaps there are better ways to utilize built-in methods in Brian to make these calculations or approximations.

import pandas as pd
import numpy as np
import brian2 as b

# Load the SWC file
swc_file_path = filtered_swc_file_path # This is a file that doesn't have axonal components
swc_data = pd.read_csv(swc_file_path, delim_whitespace=True, comment='#',
                       names=['n', 'type', 'x', 'y', 'z', 'radius', 'parent'])

# Identify the soma segment
soma_segment = swc_data[swc_data['parent'] == -1]
soma_coords = soma_segment[['x', 'y', 'z']].values[0]
soma_radius = soma_segment['radius'].values[0]

# Function to compute the Euclidean distance considering y as height
def euclidean_distance(coord1, coord2):
    return np.sqrt((coord1[0] - coord2[0])**2 + (coord1[1] - coord2[1])**2 + (coord1[2] - coord2[2])**2)

# Compute the distances of each segment from the soma
swc_data['distance_from_soma'] = swc_data.apply(
    lambda row: euclidean_distance(soma_coords, row[['x', 'y', 'z']].values), axis=1)

# Define the distance thresholds
proximal_threshold = 40 * b.um # Contains soma and dendrites close to soma
middle_threshold = 240 * b.um # Contains most basal dendrites and proximal part of apical dendrites
distal_threshold = 800 * b.um # Contains apical dendrites

# Initialize surface area counters
total_surface_area = 0.0
proximal_surface_area = 0.0
middle_surface_area = 0.0
distal_surface_area = 0.0

# Function to calculate segment distance from the soma
def segment_distance(section):
    if section.parent is None:
        return 0 * b.um
    return np.sqrt((section.x[0] - section.parent.x[0])**2 +
                   (section.y[0] - section.parent.y[0])**2 +
                   (section.z[0] - section.parent.z[0])**2)

# Function to traverse and calculate areas
def traverse_and_calculate_area(section, current_distance):
    global total_surface_area, proximal_surface_area, middle_surface_area, distal_surface_area
    segment_area = np.sum(section.area)
    total_surface_area += segment_area

    # Calculate the new distance
    segment_dist = segment_distance(section)
    new_distance = current_distance + segment_dist

    # Classify the segment area based on the distance thresholds
    if new_distance <= proximal_threshold:
        proximal_surface_area += segment_area
    elif new_distance <= middle_threshold:
        middle_surface_area += segment_area
    elif new_distance <= distal_threshold:
        distal_surface_area += segment_area

    # Recursively process child segments
    for child in section.children:
        traverse_and_calculate_area(child, new_distance)

# Start the traversal from the root (soma) and calculate areas
traverse_and_calculate_area(morpho, 0 * b.um)

print('Total surface area:', total_surface_area)
print('Proximal surface area:', proximal_surface_area)
print('Middle surface area:', middle_surface_area)
print('Distal surface area:', distal_surface_area)

Many thanks in advance for any insight on these matters.

Best wishes,
Rares

Hi @raresdorcioman. Manipulating/traversing existing morphologies is not super-convenient in Brian at the moment, but Iā€™d like to change this. Your feedback/use cases are therefore very welcome!

Regarding your first question: this looks very reasonable to me. If you want to use the morphology for simulation afterwards, you might also consider correcting the surface area of the dendrites. I am not an expert on this, but neurons often have spines on the dendrites which are not part of the reconstruction. The surface area in these reconstructions is therefore underestimated, and usually people who do simulations increase the length of the sections by a certain factor to get a more realistic surface area.

Regarding question 2: no, I donā€™t think how this was meant. When people talk about the ā€œfirst 40Āµm of dendritesā€, they are usually considering the 1-dimensional distance along the path from the soma, and not the distance in 3D space. In a Brian morphology, this is the distance attribute. For example (no 3D positions at all used here):

>>> soma = Soma(10*um)
>>> soma.dendrite = Cylinder(5*um, length=100*um, n=10)
>>> print(soma.dendrite.distance)
[ 5. 15. 25. 35. 45. 55. 65. 75. 85. 95.] um

You can see that the midpoint of the first segment of the dendrite is at 5Āµm. By convention, this is the distance from the soma, but you could add half of the somaā€™s diameter to get the distance from the center of the soma, for example.

When you use a morphology to create a SpatialNeuron, you can actually access these things more conveniently across the full morphology. For example, take this simple morphology with two dendrites:

soma = Soma(10*um)
soma.dend1 = Cylinder(5*um, length=100*um, n=10)
soma.dend2 = Cylinder(5*um, length=100*um, n=10)
neuron = SpatialNeuron(soma, 'Im = 0*nA/cm**2 : amp/meter**2')

To get the indices of all compartments at less than 45Āµm from the soma, you can use Brianā€™s string-based syntax:

>>> print(neuron.indices['distance < 45*um'])
[ 0.  1.  2.  3.  4. 11. 12. 13. 14.]

(I just noted that there is a bug ā€“ these should be integers and not floatsā€¦ You can use .astype(int) to convert them)

You can also use this to directly ask for the area, for example:

>>> print(neuron.area['distance < 45*um'])
[314.15926536 157.07963268 157.07963268 157.07963268 157.07963268
 157.07963268 157.07963268 157.07963268 157.07963268] um^2
>>> print(sum(neuron.area['distance < 45*um']))
1570.79632679 um^2

All this has a few rough edges, but I hope it gets you going. Let us know if you run into any issues!

Hi @mstimberg,

Thanks a lot for your prompt response and for the guidance on how to use the SpatialNeuron to simplify these calculations. And good point on the dendritic spines, their contribution might be quite significant in some cases, so Iā€™ll make sure to investigate this.

This looks much more manageable now that Iā€™m aware of this built-in method to compute the areas I am interested in. Iā€™m trying to get accustomed to working with such neurons, but itā€™s great to know about this, so thank you for taking the time.

Best wishes,
Rares

1 Like